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CHAPTER  I 


INTRODUCTION 

Coastal  regions  are  a  limited  resource  that  must  be  shared  among 
sometimes  conflicting  environmental,  commercial,  industrial,  residen¬ 
tial,  and  recreational  needs.  These  competing  interests  make  it  diffi¬ 
cult  and  expensive  to  design,  build,  and  repair  coastal  projects.  In 
order  to  more  efficiently  utilize  our  coastal  areas,  it  is  imperative  to 
develop  the  best  possible  methods  of  engineering  analysis.  Accurate  and 
efficient  shallow  water  wave  transformation  models  are  one  category  of 
tool  needed  for  coastal  engineering  analysis.  These  computer  models 
predict  changes  in  the  height  and  the  direction  of  waves  propagating 
through  water  of  varying  depth.  The  purpose  of  this  report  is  to  de¬ 
velop  a  numerical  model  of  the  transformation  of  finite-amplitude  waves, 
which  will  fill  a  gap  in  the  engineering  analysis  of  shallow  water 
waves . 

The  unique  features  of  this  numerical  model  include  the  following 
items:  a)  A  second  approximation  to  cnoidal  wave  theory,  a  finite- 
amplitude  theory  valid  for  intermediate  and  shallow  water,  is  used. 

(The  first  approximation  to  cnoidal  wave  theory  can  also  be  used.) 
b)  Expressions  for  the  efficient  calculation  of  cnoidal  waves  (Isobe, 
1985)  are  rederived  and  applied  in  the  numerical  model,  c)  Wave  heights 
and  angles  of  propagation  resulting  from  the  shoaling  and  refraction  of 
finite-amplitude  waves  are  calculated  directly  on  a  numerical  grid  (as 
opposed  to  the  more  cumbersome  classical  wave  ray  method). 

The  foundation  of  any  wave  transformation  model  is  the  wave  theory 
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upon  which  it  is  based.  The  two  main  families  of  perturbation  wave 
theory  are  Stokes  theory  (Stokes,  1847)  and  cnoidal  theory  (Korteweg  and 
de  Vries,  1895).  Both  are  based  on  perturbation  solutions  to  the  non¬ 
linear  equations  of  the  water-wave  boundary  value  problem.  (A  discus¬ 
sion  of  water-wave  theory  is  given  in  Chapter  II.)  Stokes  theory,  valid 
for  intermediate  and  deep  water,  is  based  on  trigonometric  functions. 
Cnoidal  theory,  valid  for  intermediate  and  shallow  water,  is  based  on 
the  Jacobian  elliptic  functions. 

The  first  approximation  to  Stokes  theory,  also  called  Airy, 
linear,  or  small-amplitude  theory,  is  the  basis  for  almost  all  wave 
transformation  models  presently  in  engineering  use.  Small-amplitude 
theory  performs  well  in  a  wide  variety  of  engineering  applications. 
However,  the  infinitesimal  wave  height  assumption  used  in  the  derivation 
of  small -amplitude  wave  theory  becomes  increasing  invalid  in  shallow 
water.  In  shallow  water,  where  the  wave  height  becomes  large  and  can 
significantly  affect  wave  properties,  cnoidal  theory  offers  an  improve¬ 
ment  over  small -amplitude  theory,  since  cnoidal  wave  theory,  even  at  the 
first  approximation,  describes  waves  of  finite  height. 

For  the  numerical  model  developed  in  this  report,  a  second  approx¬ 
imation  to  cnoidal  wave  theory  (Isobe  and  Kraus,  1983b)  is  used. 

(Chapter  II  contains  an  outline  of  the  derivation  of  this  theory.) 

The  second  approximation  was  included  to  increase  the  accuracy  of  the 
calculation  of  shallow  water  waves.  In  particular,  the  second  approx¬ 
imation  more  realistically  predicts  the  horizontal  water  particle  veloc¬ 
ity.  (The  first  approximation  to  cnoidal  theory  calculates  velocities 
which  do  not  vary  with  depth.)  Including  the  second  approximation 


produces  more  realistic  predictions  of  the  depth  dependence  of  the 
horizontal  water  particle  velocities. 

Although  cnoidal  wave  theory  has  been  available  for  over  90  years, 
its  use  in  engineering  analysis  has  been  minimal.  This  i3  because  of 
the  complexity  and  unfamiliarity  of  the  mathematics  of  cnoidal  wave 
theory,  and  also  the  difficulty  encountered  in  its  computer  calculation. 
A  goal  of  this  report  is  to  promote  cnoidal  wave  theory  as  a  viable 
engineering  design  tool.  This  will  be  done  by  providing  a  discussion  of 
the  mathematical  basis  of  cnoidal  theory  (Chapter  III),  and  by  deriving 
methods  (Appendix  A)  for  the  efficient  calculation  of  cnoidal  wave 
theory,  based  on  the  work  of  Isobe  (1985). 

Many  of  the  small -amplitude  wave  models  and  all  of  the  known 
finite-amplitude  wave  models  calculate  refractive  effects  using  ray 
theory  (Munk  and  Authur,  1952).  However,  ray  theory  is  difficult  to 
use,  and  complicated  interpolation  is  needed  in  order  to  specify  model 
results  on  a  numerical  grid.  A  method  of  calculating  results  directly  at 
grid  locations  (introduced  for  small -amplitude  theory  by  Noda  et  al., 
1974)  is  adapted  here. 

Shoaling  and  refraction  are  the  traditional  first  steps  in  the 
calculation  of  wave  transformation.  (Wave  transformation  is  discussed 
in  Chapter  IV.)  The  model  presented  in  this  report  will  calculate  the 
effects  of  shoaling  and  refraction  for  the  region  from  just  seaward  of 
the  surf  zone  out  to  water  of  intermediate  depth.  Before  the  computer 
became  an  indispensable  engineering  analysis  tool,  shoaling  and 
refraction  were  predicted  using  hand  and  graphical  calculations. 

Usually,  simplifications  were  made  to  the  bathymetry  by  assuming  a 


locally  plane  bottom  with  straight  and  parallel  depth  contours.  To  use 
a  more  realistic  bathymetry  tedius  computations  were  often  necessary. 
The  Shore  Protection  Manual  (SPM)  (SPM,  1984)  provides  tables  for 
shoaling  calculations  and  a  template  which  speeds  the  hand  calculation 
of  refraction.  With  the  advent  of  the  computer,  numerous  numerical 
transformation  models  for  irregular  bathymetry  were  developed.  The  SPh 
(SPM,  1984,  pp.  2-71)  lists  numerous  sources.  These  models  are  all 
based  on  small-amplitude  wave  theory  and  most  use  ray  theory  for  the 
calculation  of  refraction. 

Recently,  there  has  been  much  research  and  development  of  wave 
transformation  numerical  models.  Theoretical  developments  for  combined 
refraction  and  diffraction  for  small -amplitude  waves  were  pioneered  by 
Berkhoff  (1972)  and  Radder  (1979).  Ebersole  (1985)  and  Lozano  and  Liu 
(1980)  are  among  those  who  have  presented  models  using  combined  wave 
refraction  and  diffraction.  These  models  improve  accuracy  in  wave 
height  and  direction  calculations,  especially  for  simulations  over 
complex  bathymetry  where  ignoring  the  effect  of  diffraction  can  result 
in  the  calculation  of  inaccurate  wave  heights  or  the  "blow  up"  of  the 
computer  run  due  to  the  convergence  of  wave  rays  (caustics).  Models 
which  incorporate  wave-current  interaction  (Liu,  1985)  and  bottom 
friction  (Liu  and  Tsay,  1985)  with  combined  refraction  and  diffraction 
have  further  advanced  the  state  of  the  art  of  the  modeling  of  shallow 
water  wave  transformation. 

The  modeling  of  waves  using  a  spectral  approach  also  has  been  an 
area  of  active  research.  Directional  spectral  models  (see  Vincent,  1982 
for  a  review)  represent  the  random  nature  of  the  sea  by  spreading  wave 


energy  over  a  two  dimensional  array  of  frequencies  and  directions. 

These  spectral  models  have  been  used  extensively  by  the  offshore  oil 
industry. 

All  of  the  efforts  mentioned  above,  while  contributing  to  the 
ability  to  model  wave  transformation,  use  small -amplitude  wave  theory. 
As  will  be  dicussed  in  Chapter  II,  small  amplitude  wave  theory  is  based 
on  the  assumption  that  the  wave  height  is  infinitesimally  small.  This 
assumption  is  clearly  violated  for  real  sea  waves  in  the  coastal  zone. 
However,  because  of  the  increased  complexity  of  finite-amplitude  theory, 
only  a  few  finite-amplitude  numerical  models  that  calculate  wave  trans¬ 
formation  over  an  irregular  bottom  have  been  developed.  Chu  (1975)  was 
a  pioneer,  developing  a  model  which  combined  Stokes  wave  theory  at  a 
third  approximation,  cnoidal  wave  theory  at  a  first  approximation,  and 
hyperbolic  wave  theory.  Headland  and  Chu  (1985)  attempted  to  couple 
small-amplitude  theory  in  deep  and  intermediate  water  with  cnoidal 
theory  in  shallow  water.  Problems  with  connecting  the  two  components 
forced  them  to  present  separate  results  for  the  deep  and  shallow  water 
models.  A  model  for  the  shoaling  and  refraction  of  finite-amplitude 
waves,  using  Stokes  theory  at  the  third  approximation  and  incorporating 
bottom  friction,  was  presented  by  Oh  and  Grosch  (1985).  However,  the 
model  was  advocated  for  application  in  shallow  water  up  to  the  breaking 
point,  which  is  outside  the  region  of  validity  of  Stokes  theory  (Kraus, 
Cialone,  and  Hardy,  1987).  No  known  finite-amplitude  transformation 
model  calculates  directly  at  grid  points.  For  the  calculation  of 
refraction,  all  three  of  the  aforementioned  finite-amplitude  models 


employed  ray  theory,  which  requires  interpolation  to  specify  information 
at  grid  locations. 

The  numerical  model  presented  in  this  report  is  a  first  step  in 
the  development  of  a  nearshore  finite-amplitude  wave  transformation 
model.  This  model  will  couple  a  model  employing  Stokes  theory  at  the 
third  approximation  (Cialone  and  Kraus,  1987)  for  use  in  deeper  water 
with  the  cnoidal  model  described  in  this  report  for  use  in  shallower 
water.  The  final  form  of  this  combined  model  will  include  additional 
wave  transformations  and  may  employ  more  sophisticated  numerical 
techniques  such  as  a  boundary  fitted  numerical  grid.  This  report  will 
concentrate  on  the  efficient  and  accurate  computer  calculation  of  both 
cnoidal  wave  theory  up  to  a  second  approximation  and  the  basic  wave 
transformation  mechanisms,  shoaling  and  refraction.  The  numerical 
schemes  presented  at  this  first  step  will  emphasize  simplicity  and 
economy  of  calculation. 

The  remainder  of  this  report  is  structured  as  follows.  Chapter  II 
contains  a  brief  discussion  of  water-wave  theories  used  in  engineering 
analysis,  and  an  outline  of  a  derivation  of  second-order  cnoidal  wave 
theory.  In  Chapter  III,  the  role  of  Jacobian  elliptic  functions  in 
cnoidal  wave  theory  is  described,  and  expressions  for  their  efficient 
calculation  are  presented.  The  derivations  of  the  governing  equations 
used  in  the  model  for  the  calculation  of  shoaling  and  refraction  are 
included  in  Chapter  IV.  Chapter  V  discusses  the  numerical  model. 

Results  from  numerical  simulations  are  presented  and  discussed  in 
Part  VI.  And,  Chapter  VII  contains  the  conclusions  of  this  study,  and 
some  suggestions  for  future  work.  Appendix  A  contains  a  derivation  of 


the  expressions  for  the  efficient  calculation  of  cnoidal  wave  theory  and 
a  short  table  of  properties  of  Jacobian  elliptic  functions.  Derivations 
of  equations  for  wave  energy,  energy  flux,  and  velocity  of  energy  propa¬ 
gation  for  cnoidal  waves  at  a  second  approximation  are  presented  in 
Appendix  B.  Appendix  C  is  a  computer  listing  of  the  numerical  model. 
Appendix  D  contains  plots  from  the  results  of  the  refraction  of  second- 
order  cnoidal  waves.  A  list  of  symbols  used  throughout  this  report  is 
contained  in  Appendix  E. 


CHAPTER  II 


SURVEY  OF  WATER  WAVE  THEORY 

The  fundamental  element  of  any  wave  transformation  model  is  the 
wave  theory  upon  which  the  calculations  are  based.  Before  outlining  a 
derivation  of  the  second-order  cnoidal  wave  theory  used  in  this  report, 
the  basic  water-wave  boundary  value  problem  will  be  presented  and  some 
of  the  major  water-wave  theories  currently  in  engineering  use  will  be 
briefly  discussed.  The  motivation  for  presenting  this  brief  survey  of 
wave  theory  is  to  delineate  the  relationship  of  cnoidal  theory  to  the 
commonly  applied  engineering  wave  theories.  For  a  more  detailed 
presentation  of  basic  wave  theory,  Ippen  (1966,  Chapters  1  and  2), 
Sarpkaya  and  Isaacson  (1981,  Chapter  4),  and  Dean  and  Dalrymple  (1984, 
Chapter  3)  are  recommended. 

2 . 1  The  Boundary  Value  Problem  for  Water  Waves 

The  various  water-wave  theories  are,  in  essence,  composed  of 
different  sets  of  assumptions  and  techniques  which  attempt  to  describe 
the  dynamics  and  kinematics  of  real  waves.  Real  water  waves  have 
several  properties  which  make  their  mathematical  specification  very 
difficult.  Water  is  a  viscous  fluid  and  real  waves  travel  over  sea 
bottoms  or  lake  beds  which  are  in  general  neither  impermeable  nor  flat. 
Because  of  these  difficulties  most  of  the  wave  theories  used  for 
engineering  purposes  begin  with  a  series  of  identical  assumptions.  The 
problem  is  formulated  as  a  boundary  value  problem  in  which  a  governing 
differential  equation  is  solved  over  the  region  of  interest  and 
satisfies  the  appropriate  boundary  conditions. 


Figure  2.1  is  the  definition  sketch  of  the  boundary  value  problem, 
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where  the  following  symbols*  are  used: 

x,z  -  Coordinate  axis  with  the  origin  at  the  still-water  level,  z 

positive  upward,  and  the  wave  moving  in  the  positive  x  direction. 

D  -  Water  depth,  defined  as  the  distance  from  the  still-water  level  to 
the  bottom. 

n  -  Distance  from  the  still-water  level  to  the  free  surface. 

H  -  Wave  height,  defined  as  the  vertical  distance  from  the  crest  to 
the  trough. 

L  -  Wavelength,  defined  as  the  distance  between  successive  crests  or 
any  two  points  af  equal  phase. 

li  -  Water  particle  velocity. 

u  -  Horizontal  component  of  the  water  particle  velocity. 

w  -  Vertical  component  of  the  water  particle  velocity. 


W 


Figure  2.1  Definition  sketch  for  the  water  wave  boundary  value  problem 


*  Symbols  used  in  this  report  are  listed  in  Appendix  E. 
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The  assumption  that  the  flow  is  two-dimensional  is  inherent  in 
Figure  2.1.  The  wave  is  moving  in  the  positive  x  -  direction,  and  there 
is  no  component  of  the  motion  in  the  y  -  direction. 

The  assumption  that  water  is  incompressible  is  common  to  all 
surface  wave  theories.  This  assumption  is  justifiable  given  that 
water's  bulk  modulus  of  elasticity  is  very  large  relative  to  the  small 
changes  in  pressures  common  to  surface  waves.  Starting  with  this 
assumption,  the  governing  equation  is  derived  from  the  continuity 
equation  for  an  incompressible  fluid  (Equation  2.1). 

V  •  IT  =  0  (2.1) 

The  second  major  assumption  is  that  the  fluid  motion  is 
irrotational.  Clearly,  the  fluid  would  be  rotational  in  the  boundary 
layers  which  exist  at  the  bottom  and  at  the  water-air  surface. 
Fortunately,  these  boundary  layers  are  thin  and  the  major  portion  of  the 
water  column  can  be  considered  as  irrotational  for  engineering  purposes. 
With  the  assumption  of  irrotationality ,  a  velocity  potential  exists  and 
the  water  particle  velocity  can  be  specified  in  terms  of  the  velocity 
potential,  4>  (Equation  2.2).  Equation  2.2  is  substituted  into  Equation 
2.1,  resulting  in  Equation  2.3.  Rewriting  Equation  2.3  results  in  the 
Laplace  equation  (Equation  2.4),  which  is  used  as  the  governing  equation 
in  the  water-wave  boundary  value  problem. 


With  the  governing  equation  established,  specification  of  boundary 
conditions  is  necessary  before  a  solution  can  be  determined.  Boundary 
conditions  are  needed  at  the  free  surface,  the  bottom,  and  lateral  sides 
of  the  region  of  interest. 

Kinematic  boundary  conditions  are  specified  at  the  bottom  and  at 
the  free  surface.  These  boundary  conditions  force  the  fluid  to  remain 
in  contact  with  the  boundaries.  If  contact  were  not  maintained,  either 
there  would  be  gaps  between  the  fluid  and  the  boundary  or  the  fluid 
would  pass  through  the  boundary  (such  as  during  wave  breaking  or  at  a 
permeable  bottom).  Thus  the  kinematic  free  surface  boundary  condition, 
KFSBC  (Equation  2.5),  and  the  kinematic  bottom  boundary  condition,  KBBC 
(Equation  2.6),  have  the  inherent  assumptions  of  non-breaking  waves  and 
an  impermeable  bottom,  respectively. 

The  KBBC  is  further  simplified  by  the  assumptions  of  a  horizontal 
bottom  that  does  not  change  with  time  (no  scour  or  deposition),  resulting 
in  Equation  2.7. 

w  =  rv  +  un  on  z  :  n  (2.5) 

w  X 

w  =  D  -  uD  on  z  =  -D  (2.6) 

w  X 

w=0  on  z  :  -D  (2.7) 

where  subscripts  denote  partial  differentiation. 
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The  water-air  surface  requires  an  additional  boundary  condition.  Rigid 
boundaries  can  support  a  pressure  differential  across  the  boundary 
|  surface.  However,  the  water-air  surface  responds  to  pressure  changes  by 

adjusting  in  elevation  to  maintain  a  zero  pressure  differential  across 
the  boundary.  (Here  the  assumption  is  made  that  surface  tension  can  be 
I  ignored.)  This  extra  boundary  condition  is  needed  because  the  location 

of  the  free  surface  is  not  known  a  priori,  buc  is  itself  part  of  the 
solution.  The  unsteady  Bernoulli  equation  evaluated  at  the  free  surface 
is  used  as  the  the  dynamic  free  surface  boundary  condition,  DFSBC 
(Equation  2.8)  with  the  pressure,  p,  assumed  to  be  zero  or  gage  pressure. 
Here  g  is  the  acceleration  of  gravity  and  p0  is  the  Bernoulli  constant. 
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(2.8) 


The  lateral  boundary  conditions  are  expressed  as  the  requirement 
that  the  wave  is  of  permanent  form  and  is  periodic  in  both  space  and 
time.  These  periodic  lateral  boundary  conditions,  PLBC,  are  given  in 
Equations  2.9  and  2.10.  Here  t  is  the  independent  variable  for  time  and 
T  is  the  period  of  the  wave  form. 


<J>(  x ,  t )  =  $(x+L,t) 


$(x,t)  =  4>(x,t+T) 


(2.9) 


(2.10) 


The  boundary  value  problem  has  been  fully  specified  by  the  above 


governing  equation  and  boundary  conditions.  The  most  common  wave 


theories  used  in  engineering  analysis  can  be  derived  by  making  further 
assumptions  and  alterations  to  this  basic  formulation.  In  conclusion, 
the  basic  assumptions  can  be  summarized  as  follows: 

1.  Water  is  incompressible. 

2.  Water  is  inviscid  and,  therefore,  is  irrotational . 

3.  The  bottom  is  horizontal  and  does  not  change  with  time. 

4.  Both  the  bottom  and  the  free  surface  are  impermeable. 

5.  Surface  tension  on  the  free  surface  is  negligible. 

6.  The  wave  is  periodic  in  space  and  time,  i.e.,  is  of  permanent  form. 

The  set  of  equations  specifying  the  boundary  value  problem  is 
summarized  as 

V2*  =0  (2.11) 

w  s  nfc  ♦  unx  on  z  =  n  (2.12) 

+  5  (u2  +  w2)  +  gZ  =  —  on  z  =  n  (2.13) 

L  d  P 

w  =  0  on  z  =  -D  (2.14) 

4>(x,t)  =  <fr(x+L,t)  (2.15) 

4>(x,t)  =  $(x,t+T)  (2.16) 


In  the  next  four  sections  the  additional  assumptions  and  modifications 
necessary  for  the  development  of  small -amplitude,  Stokes,  cnoidal, 


solitary,  and  numerical  wave  theories  will  be  presented. 


2.2  Small-Amplitude  Wave  Theory 

A  closed  form  solution  of  the  wave  boundary  value  problem  speci¬ 
fied  by  Equations  2.11  to  2.16  is  extremely  difficult  to  find.  This 
difficulty  results  from  the  two  free  surface  boundary  conditions  which 

p  p 

contain  the  nonlinear  terms,  unx  and  1/2(u  ♦  wc).  These  nonlinear 

terms  can  be  shown  to  be  much  smaller  than  the  other  terms  in  these  two 
equations  using  the  assumption  that  the  wave  height  i3  infinitesimal. 
This  process  demands  that  both  H/D  <<  1  and  H/L  <<  1.  Neglecting  the 
nonlinear  terms,  the  KFSBC  and  DFSBC  can  be  written  as  Equations  2.17 
and  2.18,  respectively. 

w=nfc  onz=0  (2.17) 

P« 

+  gZ  =  —  on  z  =  0  (2.18) 

t  p 

The  system  of  linear  equations  (2.11,  2.17,  2.18,  2.14,  2.15, 

2.16)  can  be  solved  using  the  technique  of  separation  of  variables.  The 
resulting  equation  for  n  is 

n  =  |  cos  2n  (i  ♦  ^  (2.19) 

Small -amplitude  theory  is  also  called  either  linear  wave  theory  because 
of  the  linearizing  assumptions,  or  Airy  wave  theory  after  the  man  who 
first  presented  it  (Airy,  1845).  Small -amplitude  wave  theory  has  been 


and  continues  to  be  the  predominant  water  wave  theory.  This  is  because 
small -amplitude  theory  is  easy  to  calculate,  and  is  sufficiently 
accurate  for  many  engineering  and  scientific  applications. 

2.3  Finite-Amplitude  Wave  Theory 

Although  small -amplitude  wave  theory  performs  adequately  for  many 
engineering  applications,  certain  important  phenomena  are  inherently 
nonlinear.  For  example,  the  last  two  decades  have  brought  a  boom  in 
offshore  engineering,  particularly  in  the  oil  exploration  industry. 

Here  the  emphasis  is  placed  on  large  wave  heights,  which  negates  a  basic 
assumption  upon  which  linear  theory  is  based.  Specifically,  H/L  is  not 
very  much  smaller  than  one.  In  relatively  shallow  water,  where  H/D  is 
not  very  much  smaller  than  one,  finite  wave  height  affects  both  wave 
refraction  and  wave  shoaling,  and  thus  the  angle  of  wave  attack  on  a 
beach  and  the  location  of  the  breaker  line.  For  these  and  any  other 
problems  for  which  finite  wave  height  is  important,  the  boundary  value 
problem  must  be  solved  using  the  set  of  nonline  r  equations  (Equations 
2.11  to  2.16). 

In  the  treatment  of  finite-amplitude  waves,  specification  of  a 
moving  coordinate  system  usually  precedes  the  solution  of  the  nonlinear 
equations.  Use  of  a  coordinate  system  with  a  motion  equal  to  the  wave 
celerity  removes  the  unsteady  terms  from  the  equations  and  greatly 
simplifies  the  solution  procedure.  However,  since  the  final  goal  is  a 
solution  in  a  fixed  coordinate  system,  the  wave  celerity  must  be 
specified  externally  of  the  solution  of  the  nonlinear  problem,  i.e.,  it 
turns  out  that  the  wave  celerity  is  not  unique.  Specification  of  the 


* 

» 


I 

I 


I 


celerity  is  usually  done  through  a  physical  consideration.  Stokes 
(1847)  offered  the  two  most  commonly  used  definitions,  called  the  first 
and  second  definitions  of  wave  celerity.  The  first  definition  states 
that  the  wave  celerity  is  specified  so  that  the  horizontal  velocity 
averaged  over  a  wave  period  is  zero  if  evaluated  at  a  point  which  is 
always  in  the  fluid.  The  first  definition  of  wave  celerity  can  be 
expressed  as 


0 


u  dx 


(2.20) 


The  second  definition  requires  that  the  celerity  be  specified  so 
that  the  mass  flux  is  zero  when  averaged  over  a  wave  period.  The  second 
definition  of  wave  celerity  can  be  expressed  as 

0-1/  /nudz  dx  (2.21) 

0  -D 


In  this  report,  the  second  definition  of  wave  celerity  is  used.  It  is 
often  stated  in  the  literature  that  finite-amplitude  waves  cause  a  mass 
flux  in  the  direction  of  wave  propagation.  While  this  may  be  true  at 
some  locations  in  the  nearshore  wave  field,  it  is  not  true  for  an 
average  over  the  wave  field  in  the  steady  state  condition.  For  the 
steady  state  condition,  it  does  not  seem  proper  to  account  for  the 
forward  movement  of  mass  without  accounting  for  the  return  current  which 
this  forward  flux  must  cause.  Imagine  the  consequences  for  coastal 
regions  if  waves  always  caused  a  net  forward  mass  flux!  Often  the 
second  definition  of  celerity  is  proposed  for  matching  theoretical 
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results  with  wave  flume  experimental  data.  The  closed  system  of  the 
wave  flume  demands  that  a  no-flux  condition  exist  after  the  steady  state 
is  reached.  Granted  that  the  natural  environment  can  be  far  more 
complex  than  a  wave  flume;  however,  a  steady  state  condition  in  nature 
requires  the  same  balance  of  inflows  and  outflows  as  does  the  wave 
flume.  Therefore,  the  zero  mass  flux  of  the  second  definition  would 
seem  to  give  the  best  overall  representation.  The  second  definition  of 
celerity  has  another,  more  practical  benefit.  For  the  derivation  of  an 
expression  to  calculate  energy  flux  (Appendix  B),  the  second  derivation 
allows  a  substantial  simplification  (see  Equations  B.9  and  B.10).  To 
keep  this  matter  in  perspective,  it  is  important  to  note  that  the  choice 
of  one  of  these  two  definitions  is  a  relatively  academic  decision.  The 
choice  does  not  effect  the  first  two  approximations  to  the  wave  celerity 
and  only  results  in  a  small  difference  at  higher  order  approximations. 

The  perturbation  method  is  used  to  solve  the  set  of  nonlinear 
equations  (2.11  to  2.16).  In  this  technique,  the  dependent  variables  in 
the  boundary  value  problem  are  assumed  to  be  functions  of  an  auxiliary 
parameter,  6  ,  and  are  expanded  in  a  power  series  of  a  perturbation 
parameter,  e  ,  which  is  assumed  to  be  small.  For  example,  Equation  2.22 
is  a  power  series  representation  for  n  . 

n  =  en^  +  e2rt2  +  +  ...  (2.22) 

Waves  of  permanent  form  can  be  described  by  three  independent 
variables,  H,  L,  and  D.  Two  independent  nondimens ional  quantities  can 
be  formed  from  these  three  variables.  The  most  commonly  used  are  H/L 


and  H/D.  One  of  these  quantities  is  usually  chosen  as  e  .  For  a 
derivation  of  Stokes  wave  theory,  H/L  is  often  chosen  as  the 
perturbation  parameter;  whereas,  for  a  derivation  of  cnoidal  wave 
theory,  H/D  is  often  chosen. 

The  power  series  representations  of  the  variables  are  inserted 
into  the  system  of  equations.  Then  the  terms  with  an  equal  power  of  c 
in  each  equation  are  grouped  together  and  solved  separately.  For 
example,  the  equations  of  first-order  involve  only  those  terms  multi¬ 
plied  by  e'  ,  the  equations  of  second-order  involve  only  those  terms 
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multiplied  by  e  ,  and  the  equations  of  nth-order  involve  only  those 
terms  multiplied  by  e11  .  The  terminology,  the  nth-order  solution  or  the 
nth-order  theory,  implies  the  solution  of  the  equations  up  to  and 
including  nth-order.  The  solution  of  the  equations  of  each  successive 
order  adds  a  correction  to  the  (n-1 )th-order  solution.  Unfortunately, 
the  higher  order  solutions  are  gained  at  considerable  algebraic  expense. 

2.4  Stokes  Wave  Theory 

Sir  George  Stokes  originated  the  family  of  wave  theories  bearing 
his  name  (Stokes,  1847)  when  he  developed  a  second-order  theory.  This 
derivation  showed  that  the  first  approximation  was  identical  to  Airy 
wave  theory;  consequently,  first-order  Stokes  theory  and  small -amplitude 
theory  are  identical.  There  have  been  several  derivations  of  higher 
order  Stokes  wave  theories.  Skjelbreia  and  Hendrickson  (I960)  obtained 
a  fifth-order  solution.  Their  derivation  used  the  same  expansion 
parameter  as  Stokes,  ak,  where  a  =  H/2  is  the  wave  amplitude  at  first 
order  and  the  wave  number,  k  =  2ir/L  .  (Nishimura,  Isobe  and  Horikawa 


(1977)  found  a  minor  error  in  this  solution.)  Derivations  using  e 
=  kH/2,  include  Isobe  and  Kraus  (1983a)  to  third-order  and  Fenton  (1985) 
to  fifth-order. 


Stokes  wave  theories  are  not  valid  over  the  whole  range  of 
relative  depths,  D/L.  Isobe,  Nishimura  and  Horikawa  (1982)  showed  that 
their  perturbation  solution  of  Stokes  theory  is  valid  only  if  the  Ursell 
parameter  (Equation  2.23)  has  a  value  less  than  25.  (For  other 
perturbation  solutions  this  limit  may  vary  slightly.)  This  implies 
validity  only  in  intermediate  to  deep  water.  If  Stokes  theory  is 
applied  outside  its  range  of  validity,  abnormalities  such  as  secondary 
peaks  of  the  wave  free  surface  profile  occur  soon  after  the  limiting 
value  of  U  is  exceeded.  Naturally,  the  numerical  values  of  certain 
physical  quantities  will  also  be  inaccurate  outside  the  range  of 
validity. 


U  = 


(2.23) 


H  3 

Goda  (1983)  presented  an  empirical  parameter,  nr—  coth3  k_.D  , 

LSA 

where  the  subscript  SA  denotes  a  quantity  calculated  from  small- 
amplitude  theory.  This  parameter  smoothly  connects  a  multiple  of  the 
Ursell  parameter,  which  is  used  as  a  measure  of  nonlinearity  in  shallow 
water  with  the  deepwater  wave  steepness,  H/L,  which  serves  as  a  measure 
of  nonlinearity  in  deep  water.  Goda  found  that  Stokes  theory  was  valid 
below  n  »  0.2  .  This  parameter  has  not  been  thoroughly  studied,  but 
further  research  might  prove  it  to  be  a  valuable  descriptor  of  finite 
amplitude  waves. 
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2.5  Cnoidal  Wave  Theory 


The  Dutch  mathematicians,  Korteweg  and  de  Vries  (1895),  during  an 
investigation  of  the  change  in  form  of  propagating  long  waves,  developed 
a  type  of  long  wave  surface  profile  which  can  be  specified  in  terms  of 
the  Jacobian  elliptic  function,  cn.  Using  an  analogy  with  waves 
described  by  sinusoidal  functions,  they  coined  the  term  "cnoidal"  for 
this  wave  type. 

Keulegan  and  Patterson  (1940),  as  part  of  treatise  on  irrotational 
waves,  developed  a  well-known  first-order  cnoidal  solution.  Friedrichs 
(1948)  introduced  a  systematic  method  for  deriving  shallow  water  wave 
theories  by  using  a  coordinate  stretching  nondimensionalization.  This 
technique  is  the  basis  for  many  of  the  succeeding  derivations  of  cnoidal 
theory.  The  most  commonly  cited  derivations  of  cnoidal  theory  include 
Keller  (1948),  Laitone  (I960),  and  Chappelear  (1962).  Svendsen  (1974), 
in  a  report  which  gives  a  thorough  discussion  of  first-order  cnoidal 
wave  theory,  provides  tabulated  values  of  the  parameters  necessary  for 
the  engineering  calculation  of  first-order  cnoidal  waves. 

The  region  of  applicability  of  cnoidal  theory  is  in  intermediate 
to  shallow  water,  up  to  but  perhaps  not  including  the  wave  breaking 
point.  Svendsen  (1974)  derived  a  deepwater  limit  for  cnoidal  theory. 
This  limit  varies  between  D/LQ  =  0.12  for  waves  of  infinitely  small  wave 
height  to  D/LQ  =  0. 14  for  waves  that  are  steep  enough  to  break  at  the 
deepwater  limit.  LQ  is  the  small-amplitude  theory  deepwater  wavelength, 
Lq  =  gT  /(2ir)  .  Isobe,  Nishimura  and  Horikawa  (1982)  investigated  the 
error  in  the  DFSBC  for  cnoidal  solutions  of  several  different  orders. 
Their  findings  indicate  a  deepwater  limit  close  to  D/L  =  0.2.  (Note 


that  the  actual  wavelength  is  used.)  These  values  of  relative  water 
depth  indicate  that  Stokes  and  cnoidal  theories  have  an  overlapping 
range  of  Ursell  numbers  in  which  both  theories  are  valid,  approximately 
25  >  U  >  10.  Cnoidal  theory  is  valid  at  higher  Ursell  numbers  and  Stokes 
theory  is  valid  at  lower  Ursell  numbers.  Because  of  the  complexity  of 
cnoidal  theory,  Isobe,  Nishimura,  and  Horikawa  (1982)  suggest  that 
Stokes  theory  be  used  in  the  region  of  overlap  in  validity  with  cnoidal 
theory.  However,  calculation  techniques  introduced  by  Isobe  (1985) 
increase  the  ease  of  cnoidal  wave  calculations  and  diminish  the 
difference  in  the  difficulty  of  calculation  between  the  two  theories. 
These  expressions  are  introduced  in  Chapter  III,  derived  in  Appendix  A, 
and  used  in  the  numerical  model  developed  in  this  thesis. 

2.6  Solitary  Wave  Theory 

Scott  Russell,  while  serving  as  an  engineering  consultant  to  an 
English  river  barge  company,  first  studied  solitary  waves  in  the  1830's. 
His  historic  report  (Russell,  1844)  describing  these  waves  of 
translation  served  as  a  catalyst  which  helped  stimulate  much  of  the 
early  theoretical  work  in  water  wave  theory.  Boussinesq  (1872), 

Rayleigh  (1876),  and  McCowan  (1891)  are  the  commonly  cited  classical 
references  for  the  mathematical  development  of  solitary  wave  theory. 

The  solitary  wave  is  related  to  cnoidal  waves,  as  the  same 
assumptions  and  equations  are  used  to  specify  the  boundary  value 
problem.  The  solitary  wave  results  as  a  limiting  case  of  the  cnoidal 
wave  as  the  wavelength,  and  therefore,  the  wave  period  approaches 
infinity.  For  a  solitary  wave,  the  wave  form  exists  totally  above  the 
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still  water  surface,  and  the  water  particle  velocities  have  only  a 
forward  motion.  The  periodic  motion  which  characterizes  the  other  wave 
theories  is  replaced  by  a  wave  of  translation  in  solitary  wave  theory. 
Since  waves  in  the  surf  zone  often  take  on  the  appearance  of  solitary 
waves,  these  waves  have  been  suggested  for  engineering  use  in  very 
shallow  water  (Munk,  1949). 

2.7  Numerical  Wave  Theories 

Increasing  the  order  of  a  perturbation-type  analytical  wave  theory 
is  done  only  with  great  effort.  Also,  two  analytical  theories  are 
needed  to  cover  the  full  range  of  water  depths.  Therefore,  methods  were 
sought  which  would  accurately  solve  the  boundary  value  problem  without 
the  need  for  the  derivation  of  a  high-order  analytical  theory.  Dean 
(1965)  presented  such  a  method  which  he  called  the  Stream  Function  Wave 
Theory.  The  KFSBC  is  automatically  satisfied  if  the  stream  function  is 
used  instead  of  the  velocity  potential.  The  only  boundary  condition  not 
exactly  satisfied  is  the  DFSBC.  (The  use  of  the  stream  function  is  not 
unique  to  this  theory,  but  is  also  commonly  used  in  both  Stokes  and 
cnoidal  theory  derivations.)  The  coefficients  of  a  Stokes-like  expansion 
are  determined  for  a  given  H,  T,  and  D,  so  that  the  least-square  error 
of  the  DFSBC  is  reduced  to  a  very  small  value.  The  theory  has  been  shown 
to  very  accurately  match  experimental  data  (Dean,  1974). 

Cokelet  (1977)  presented  a  computer-generated  solution  of  a  very 
high  order  Stokes-like  expansion.  He  used  a  specially  developed 
perturbation  parameter  and  Pade  approximates,  a  mathematical  tool  for 
summing  series.  This  solution  is  considered  an  exact  solution  to  the 


nonlinear  water  wave  boundary  value  problem  for  waves  of  permanent 
form.  Results  obtained  using  Cokelet's  theory  are  often  used  for 
comparison  with  the  results  from  other  theories.  Nishimura,  Isobe  and 
Horikawa  (1977)  provide  very  high  order  solutions  to  both  cnoidal  and 
Stokes  waves  using  computer  solved  recurrence  formulae.  They  used  these 
solutions  to  study  the  accuracy  with  which  the  higher  order  theories 
matched  the  free  surface  boundary  conditions. 

These  and  other  computer-dependent  solution  techniques  are 
important  for  the  theoretical  study  of  such  problems  as  determining 
regions  of  validity,  and  for  engineering  studies  where  the  problem  is 
highly  nonlinear,  such  as  calculating  the  forces  from  the  highest 
possible  wave.  However,  a  major  disadvantage  results  from  these 
computer  generated  solutions.  There  are  no  explicit  formulae  for 
engineering  quantities,  so  every  combination  of  H,  T,  and  D  requires  a 
separate  solution.  Tabulated  values  have  been  provided  of  the  important 
engineering  quantities  for  various  values  of  H,  T,  and  D  (Dean,  1 974 ) . 
However,  use  of  either  the  tables  or  the  numerical  techniques  would  be 
difficult  and  expensive  in  a  wave  transformation  model  where  thousands 
of  different  combinations  of  the  basic  parameters  would  result  for  each 
simulation. 

2.8  Brief  Outline  of  a  Derivation  Second-Order  Cnoidal  Wave  Theory 

The  finite-amplitude  wave  transformation  described  in  this  report 
is  intended  for  use  in  intermediate  and  shallow  water,  seaward  of  the 
breaker  line.  Stokes  wave  theory  is  not  applicable  to  this  region. 
Solitary  wave  theory  does  not  adequately  describe  the  periodic  nature  of 
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the  waves  in  the  region  outside  of  the  surf  zone.  Numerical  theories 
are  not  computationally  efficient.  Therefore,  cnoidal  theory  was 
chosen.  Second-order  theory  was  selected  to  increase  the  accuracy  with 
which  the  theory  would  predict  the  needed  engineering  quantities.  In 
particular,  second-order  theory  more  realistically  predicts  the  vertical 
variation  of  the  horizontal  water  particle  velocity,  since  first-order 
theory  predicts  a  uniform  velocity  from  the  water  surface  to  the  sea 
bottom. 

The  derivation  of  second-order  cnoidal  theory  (Isobe  and  Kraus, 
1983b)  was  chosen  for  two  reasons.  First,  this  derivation  is  rigorous, 
systematic,  and  easy  to  use.  Second,  a  companion  derivation  for  third- 
order  Stokes  theory  (Isobe  and  Kraus,  1983a)  provides  the  foundation  for 
a  numerical  model  which  will  serve  as  a  companion  refraction  model, 
providing  seaward  boundary  conditions  for  the  model  developed  in  this 
report. 

It  is  not  the  purpose  of  this  report  to  fully  document  the  second- 
order  cnoidal  solution.  A  brief  outline  of  the  derivation  will  be 
presented  in  order  to  provide  some  background  for  the  reader. 

Using  the  relationship  between  the  velocity  potential  and  stream 
function  (Equation  2.24),  and  assuming  a  coordinate  system  which  moves 
with  the  speed  of  the  wave  form,  the  governing  equation  and  boundary 
conditions  for  the  nonlinear  wave  boundary  value  problem  (Equations  2.11 
to  2.16)  are  replaced  by  Equations  2.25  to  2.30. 

i>z  =  4>x  =  u  ;  =  4»z  =  w  (2.24) 
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(2.27) 


ii»  =  0 


n  =  0 


on  z  =  -D  (2.28) 


(2.29) 


where  the  over  bar  denotes  the  quantity  is  averaged  over  wavelength,  and 


n(0)  -  r\(L/2)  =  H 


(2.30) 


Note  that  the  use  of  the  stream  function  simplifies  the  KFSBC  (Equation 
2.26),  as  was  mentioned  in  discussing  Dean's  Stream  Function  Wave 
Theory.  Also  the  PLBC  (Equations  2.15,  2.17)  are  replaced  by  Equations 
2.29  and  2.30.  Equations  2.25  to  2.30  can  be  nondimensionalized  using 
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Substituting  Equations  2.31  to  2.37  into  Equations  2.2 5  to  2.30  results 
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on  Z  =  N 
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f  =  0 


on  Z  =  -1 


(2.41) 


N  =  0 


(2.42) 


N(0)  -  N(1/2)  =  g 


(2.43) 


The  term,  D2/L  ,  which  appears  in  Equations  2.38  and  2.40,  as  will  be 


seen  more  clearly  below,  causes  a  difference  in  order  between  the  x  and 


z  derivatives.  It  is  this  feature  which  determines  that  these  equations 


will  produce  a  shallow  water  wave  solution. 


After  the  equations  are  readied,  the  next  step  in  setting  up  the 
perturbation  method  of  solution  is  to  select  an  ordering  parameter.  For 
this  derivation  the  ordering  or  perturbation  parameter  is 


e 


(2.44) 


Next,  all  the  variables  in  the  Equations  2.38  to  2.43  are  expressed  as 
functions  of  the  auxliliary  parameter,  6  ,  and  expanded  as  a  power 
series  about  e  .  In  cnoidal  theory,  the  modulus  of  the  elliptic 
integral,  <  ,  is  usually  selected  as  the  auxiliary  parameter.  These 
expansions  are  given  in  Equations  2.45  to  2.49. 


(X,Z,K,e)  =  *0  +  ef1  +  e2¥2  +  e3^  +  ... 

(2.45) 

N(X ;k, e)  =  eN1  +  e2N2  +  e3N3  +  ... 

(2.46) 

Q(ic,e)  =  Qq  +  eQ1  +  e2Q2  +  e3Q3  +  ... 

(2.47) 

PBU,e>  =  P0  .  eP,  ♦  «2P2  *  e3P3  *  ... 

(2.48) 

(c)  1  ed1  *  c%  *  'S  * 

(2.49) 

Note  that  all  the  series  in  Equations  2.45  to  2.49  begin  with  a  zeroth- 

P  P 

order  quantity  except  the  expansions  of  N  and  D  /L  ,  which  start  with 


the  first  order.  N  is  always  smaller  than  e  ;  therefore,  N  will  not 


have  a  zeroth-order  term.  Some  other  derivations  of  cnoidal  theory, 
including  Keller  (1948),  Laitone  (I960),  and  Chappelear  (1962),  have 
used  Dc/L  as  the  perturbation  parameter.  These  two  choices  of 
perturbation  parameter,  H/D  and  D  /Lc  are  related  by  the  Ursell 
parameter  (Equation  2.23)  as  shown  in  Equation  2.50.  Since  cnoidal 
theory  is  valid  for  Ursell  numbers  greater  than  approximately  ten  (as 
was  discussed  in  Section  2.5),  Dc/L  is  smaller  than  e  ,  and  DVL 
cannot  have  a  zeroth-order  term  when  expanded  in  a  power  series 
about  e  . 


2  2 
D  /Lr 


(2.50) 


The  expanded  form  of  the  variables  are  inserted  into  the  Equations  2.38 
to  2.43.  For  example,  if  Equation  2.45  is  substituted  into  the  Equation 
2.38  the  results  are 


i^ozz  +  e’lZZ  +  e  f2ZZ  +  e  f3ZZ  +  ***^  +  ^Ed1  +  E  d2 

+  e3d3  +  ...^  x  |l-oxx  +  ef1xx  +  e  *2XX  +  ...^  =  0 


+  e  d^  +  . 


(2.51) 


If  the  wave  height  were  zero,  the  only  motion  would  be  the  uniform 
horizontal  flow  resulting  from  the  moving  coordinate  system.  Therefore, 
the  zeroth-order  contribution  to  the  stream  function  can  be  expressed  as 
Equation  2.52,  where  bQQ  is  a  constant  to  be  determined.  Equation  2.53 
follows  directly  from  Equation  2.52. 


y  =  y  =  0 
CXX  OZZ 


(2.53) 


According  to  the  perturbation  method,  terms  are  grouped  in  equal  powers 
of  e  .  The  terms  of  each  order  on  the  right  side  of  an  equation  must 
equal  the  terms  of  like  order  on  the  left  side  of  the  equation.  With 
this  grouping,  Equation  2.51  becomes 
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(2.54) 
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The  boundary  conditions  (Equations  2.39  to  2.43)  are  treated  in  a 
similar  manner.  However,  special  treatment  is  needed  for  the  two  free 
surface  boundary  conditions.  Since  the  location  of  the  free  surface  is, 
a  priori,  unknown,  variables  evaluated  at  the  free  surface  are  expanded 
using  a  Taylor  series  expanded  about  z  =  0.  For  example,  the  expansion 
for  4i(x,n)  is 


y(X,N)  =  y(X,0)  +  N  yz(X,0)  +  fzz  (X,0)  ♦  ...  (2.57) 


After  the  expanded  variables  are  inserted  into  the  governing  equation 
and  all  five  boundary  conditions  and  the  coefficients  of  each  power  of 
e  are  separated,  a  set  of  equations  for  each  order  results.  For 

p 

illustration,  the  equations  of  second-order  (coefficients  of  e  )  are 
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V2  =  0 


on  Z  =  -1 


(2.61) 


n2  =  0 


(2.62) 


N2(0)  -  N2( 1/2)  =  0 


(2.63) 


It  is  not  possible  to  obtain  the  solution  of  the  nth-order  equations 
without  using  the  equations  of  (n+1 )th-order.  The  solution  is  not 
closed  at  the  level  of  each  order,  but  instead  spirals  upward  to  include 
information  from  the  next  order.  The  solution  procedure  can  be 
summarized  as  follows.  The  equations  of  zeroth-order  will  determine  PQ 
and  Qq,  once  6qq  is  determined.  The  equations  of  first-order  determine 


b qq  so  that  a  nontrivial  solution  for  and  ^  will  exist.  Then,  if  a 


nontrivial  solution  for  n 2  and  *2  exists,  and  can  be  determined 


from  the  equations  of  second  order.  Likewise  the  equations  of  third 


order  are  used  to  determine  n2  and  y2  . 


As  was  mentioned  in  Section  2.4,  the  artifact  of  a  moving 
coordinate  system  requires  that  the  wave  celerity  be  specified  from 
considerations  other  than  the  governing  equations  and  boundary 
conditions.  The  derivation  used  in  this  report  employed  the  second 
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definition  of  Stokes.  With  the  mass  flux  equal  to  zero,  the  celerity 
becomes 

C  =  -  g  (2.64) 

With  Equation  2.64  and  after  re-dimensional izing  the  variables,  the 
solution  of  second-order  cnoidal  theory  in  a  fixed  coordinate  system 
gives 

0  ■  VsEr(co  eC1  *  z\) 


where  from  Appendix  A  (Equations  A. 76  and  A. 84) 


x 


(2.70) 

(2.71) 


and  K  (Equation  A. 2)  and  Eg  (Equation  A. 4)  are  the  complete  elliptic 
integrals  of  the  first  and  second  kind,  respectively, 


and 


e(X  - 

u>  ♦  e2(~2>  *  f  -  2lp) 

(2.72) 

A,  =  c  .  ,2  (-  |  ) 

(2.73) 

(2.74) 

=  e(X 

-  »>  ♦ 

(2.75) 

B10 

•  •  *  »2(— r^) 

(2.76) 

B20  •  c2  (-1) 

(2.77) 

(2.78) 

Bn  =  e2  (3  -  3X) 

(2.79) 
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B21  =  e2  (~  §  )  (2.80) 

e  =  £  (2.81) 

c,  .  cn2[2K(i-i)]  (2.82) 

°2  •  <="*  [2K  (c  -  f)]  <2'83> 

C3d  =  on[2i((i  -  f)]«n[a(c  '  l)]dn[2lt(c  '  f)]  (2-81,> 


and  cn,  sn,  and  dn  are  Jacobian  elliptic  functions  with  modulus  k  . 

Equation  2.69,  the  dispersion  relationship  for  second-order 
cnoidal  wave  theory,  can  be  more  conveniently  expressed  if  it  is  put  in 
terms  of  T  rather  than  L.  From  the  relation  C  =  L/T  and  using  Equation 
2.65,  Equation  2.69  can  be  rewritten  as 


16k2K2 


(2.85) 


These  results  are  compared  to  the  results  of  first-  and  second-order 
cnoidal  theories  and  linear  theory  for  H=2m,  T=8s,  and  D  =  4  m 
(U  =  81)  in  Figures  2.2  to  2.4.  As  seen  in  Fif are  2.2,  the  free-surface 
profiles  of  the  two  cnoidal  theories  are  very  similar,  and  differ  sig¬ 
nificantly  from  that  of  linear  theory.  The  asymmetry  of  the  surface 
profiles  also  appears  in  the  plots  of  the  horizontal  water  particle 
velocities  at  the  surface  and  bottom,  Figures  2.3  and  2.4, 
respectively.  (This  asymmetry  of  water  particle  velocities  is  very 
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Figure  2.2.  Free  surface  profiles  of  small -amplitude,  first-order 
cnoidal,  and  second-order  cnoidal  waves. 
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Figure  2.3.  Horizontal  water  particle  velocity  at  z  =  -D  for  small- 
amplitude,  first-order  cnoidal,  and  second-order  cnoidal  waves. 
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Figure  2.4.  Horizontal  water  particle  velocity  at  z  =  n  for  small 
amplitude,  first-order  cnoidal,  and  second-order  cnoidal  waves. 


important  for  an  accurate  description  of  sediment  transport.)  It  is 
interesting  to  note  that,  in  this  example,  although  the  surface  profiles 
of  the  first  and  second-order  cnoidal  waves  are  almost  identical,  the 
horizontal  velocities  differ  significantly.  As  in  small-amplitude 
theory  applied  to  shallow  water,  first-order  cnoidal  theory  predicts  no 
vertical  variation  in  horizontal  velocity,  whereas  second-order  cnoidal 
theory  predicts  a  greater  velocity  at  the  surface  than  at  the  bottom. 

As  can  be  seen  from  the  above  results,  the  use  of  cnoidal  theory 
depends  upon  the  capability  to  calculate  Jacobian  elliptic  functions. 
Chapter  III  will  attempt  to  provide  understanding  of  the  role  these 
functions  have  in  cnoidal  theory,  as  well  as  present  expressions  for 


their  efficient  calculation. 


CHAPTER  III 


USE  OF  ELLIPTIC  FUNCTIONS  IN  CNOIDAL  WAVE  THEORY 

3.1  An  Illustration  of  the  Role  of  the  cn  Function  in  Cnoidal  Wave 
Theory 

Most  coastal  engineers  and  coastal  scientists  have  little  or  no 
experience  with  elliptic  integrals  or  Jacobian  elliptic  functions.  This 
lack  of  familiarity  with  the  mathematical  basis  of  cnoidal  wave  theory 
is  a  hinderance  to  its  dissemination  and  use  as  an  analysis  tool.  Some 
insight  can  be  gained  into  the  role  of  Jacobian  elliptic  functions  in 
cnoidal  wave  theory  by  comparing  the  surface  profiles  of  cnoidal  and 
small -amplitude  waves.  The  following  is  not  intended  to  be  a  rigorous 
discussion  of  elliptic  functions,  but  it  will  serve  as  an  introduction 
to  the  mathematical  basis  of  cnoidal  wave  theory.  For  a  more  complete 
treatment  of  elliptic  functions,  Byrd  and  Freidman  (1954),  Milne-Thomson 
(1950),  and  Neville  (1951)  are  recommended. 

Figure  3.1  compares  the  free  surface  profiles  of  small -amplitude 
and  first-order  cnoidal  waves.  The  two  profiles  have  obvious 
differences.  The  small-amplitude  wave  is  symmetric  about  the  still- 
water  level,  whereas  the  cnoidal  wave  is  not.  The  crest  of  the  cnoidal 
wave  is  more  peaked  and  reaches  a  higher  elevation,  and  the  trough  is 
longer,  flatter,  and  at  a  higher  elevation  than  the  corresponding  parts 
of  a  small-amplitude  wave.  Assymmetry  of  the  cnoidal  wave  profile  is 
expected,  since  waves  in  the  nearshore  zone  (shallow  water)  typically 
have  sharp  peaks  and  long  shallow  troughs. 
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Figure  3.1.  Cnoidal  wave  and  small -amplitude  wave  profiles 


Although  Figure  3.1  shows  that  there  are  distinct  differences 
between  small-amplitude  and  cnoidal  waves,  an  examination  of  the 
equations  for  the  wave  profiles  will  illustrate  some  important 
similarities.  Equations  3.1  and  3.2,  the  equations  for  the  small 
amplitude  and  first-orden  cnoidal  wave  profiles,  respectively,  are 
almost  directly  analogous: 

"SA  =  I  C03[2’(c  -  !)] 

=  I  cos  0SA  (3.1) 


n  =  H  jcn2  [2K  (*  -  i)]  -  c„2  | 


■•y 

=  H  (  cn2  0  -  cn2  j 

(3.2) 

4" , 
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a.  Wave  height  is  used  in  the  amplitude  portion  of  both 
equations. 

b.  Periodic  functions,  cosine  and  cn2  ,  are  the  foundation  of 
both  equations. 

c.  The  phase  functions  are  similar.  The  only  difference  is  that 
cosine  is  periodic  with  2k  ,  whereas  cn  is  periodic  with 
2K  .  (The  quantity,  K  ,  Equation  A. 2,  is  the  complete 
elliptic  integral  of  the  first  kind.) 

d.  There  are  differences  in  structure  between  the  two  equations, 

2  p 

the  cn  ,  the  value  of  cnce  averaged  over  wavelength  (see 

Equation  A. 77),  and  the  absence  of  the  factor  of  1/2  in  the 

equation  for  the  cnoidal  wave  profile.  These  differences 

result  because  the  range  of  cosine  varies  from  -1  to  1, 

whereas  the  range  of  cn2  (as  will  be  illustrated  below) 

varies  from  0  to  1 .  Since  the  range  of  cosine  is  twice  that 

of  cn2  ,  the  "missing"  1/2  in  Equation  3.2  maintains  the  same 

? 

wave  height  for  both  profiles.  The  values  of  cn  8  are 

2 

always  positive;  therefore,  the  cn  term  is  needed  to 
position  the  cnoidal  profile  at  the  still-water  level. 

e.  The  periodic  functions,  cosine  for  the  small -amplitude  waves 
and  cn  for  cnoidal  waves,  are  the  most  fundamental 
differences  between  Equations  3.1  and  3.2.  However,  even  in 
this  most  basic  difference  there  is  a  direct  connection 
between  the  equations  for  the  two  wave  profiles.  The  cn 
function  is  defined  in  terms  of  the  cosine,  as 

cn  5  =  cos  $  (3.3; 
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The  elliptic  integral,  5  ,  is  defined  as 


/(- 


(3.4) 


2  .  2 
k  sin  6 


The  creation  of  cn  £  from  cos  $  will  be  explained  in  the 
following  discussion  and  illustration,  which  are  modified  from  Isobe  and 
Kraus  (1983b,  Appendix  B).  Each  phase,  C  ,  of  the  cn  function 
(Figure  3.2)  corresponds  to  a  value  of  the  elliptic  integral  (step  1). 
This  value  of  5  uniquely  determines  a  phase,  $  ,  of  the  cosine 
function  (step  2).  The  value  of  cos  $  (step  3)  defines  the  value  of 
cn  £  as  in  Equation  3.3  (step  4).  Note  that  cn  f;  is  periodic  with 
4K  and  has  values  ranging  from  -1  to  1.  By  squaring  cn£  ,  the 
quantity  cn  £  (the  shallow  water  wave  shape)  is  created  with  a  period 
of  2K  and  a  range  of  values  from  0  to  1 .  In  essence,  the  elliptic 
integral  is  used  as  a  mapping  function  that  "stretches"  the  symmetrical 
cosine  shape  into  the  elongated  trough  and  sharper  peak  of  the  cn*- 
shape.  Thus  it  can  be  seen  that  mathematics  of  cnoidal  wave  theory  is 
not  as  radical  a  departure  from  the  more  familiar  mathematics  of  small- 
amplitude  wave  theory  as  might  first  appear. 

An  interesting  feature  of  Jacobian  elliptic  functions  is  that  they 
are  periodic  with  a  multiple  of  K  ,  which  is  not  a  constant.  Unlike 
the  cosine  which  is  periodic  with  the  multiple  of  a  constant,  it  ,  the 
K  of  the  Jacobian  elliptic  functions  varies  with  the  modulus,  tc  . 
Figure  3.3  (adapted  from  Isobe  and  Kraus  1983b)  shows  the  effect  of 
changes  in  k  on  the  values  of  K  .  Since  k  (from  Equation  2.69) 


varies  with  the  value  of  the  Ursell  number  (Equation  2.23),  the  shape  of 
the  wave  profile  is  altered  as  the  Ursell  number  changes  value.  For 
Ursell  numbers  10  and  500,  K  is  approximately  1.88  and  9.68, 
respectively.  At  large  Ursell  numbers,  the  wave  profile  approaches  a 
solitary  wave  form.  At  small  Ursell  numbers,  the  profile  tends  more 
toward  the  shape  of  a  small-amplitude  wave.  Thus  cnoidal  wave  theory 
encompasses  the  major  wave  forms  observed  in  shallow  water. 

3.2  Efficient  Calculation  of  Elliptic  Quantities 

Methods  of  calculating  elliptic  integrals  and  Jacobian  elliptic 
functions  are  neither  as  readily  available  nor  as  easy  to  execute  as 
are  the  methods  of  calculating  the  trigonometric  functions  of  small- 
amplitude  wave  theory.  This  has  been  a  major  drawback  to  the  use  of 
cnoidal  wave  theory. 

At  an  Ursell  number  of  25  the  value  of  k  is  approximately  0.92. 
As  the  Ursell  number  increases,  k  approaches  1.0.  Subroutines  in 
computer  mathematics  libraries  and  formulae  found  in  mathematics 
handbooks  for  the  calculation  of  elliptic  functions  are  not  conveniently 
adapted  to  this  very  small  range  of  k  values.  Consequently,  use  of 
these  routines  or  formulae  can  result  in  more  expensive  computer 
calculations  and  perhaps  inaccuracies  in  the  quantities  calculated. 

Isobe  (1985)  introduced  expressions  for  the  efficient  calculation  of 
cnoidal  waves.  (Since  the  derivation  of  these  expressions  does  not 
appear  to  have  been  published  in  English,  their  development  (Isobe, 
private  communication,  1985)  is  included  in  Appendix  A.)  These  formulae 
are  given  as  power  series  of  the  complementary  nome  of  the  theta 


function,  q'  .  The  definition  of  q'  is  repeated  here  from 
Equation  A. 7  as 


K 

-1V 1 

q'  =  e  (3.5 

The  maximum  value  of  q'  applicable  to  cnoidal  theory  (U  =  10)  i 
q'  a  0.04  .  Therefore,  as  can  be  seen  in  the  following  equations,  the 
power  series  converge  very  rapidly,  and  only  the  first  two  or  three 
terms  are  needed  to  ensure  accuracy.  These  efficient  expressions  are 
the  bases  of  the  cnoidal  wave  calculations  in  the  numerical  model  which 
is  developed  in  this  thesis. 

The  results  of  the  derivation  of  second-order  cnoidal  theory 
presented  in  the  previous  chapter  show  that  expressions  are  needed  for 
the  following  seven  quantities  (Equations  3.6  to  3.12).  This 
presentation  of  equations  follows  from  Isobe  (1985). 
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where 


(T03)‘ 


p  c. 

( B' )  =  sinh  S'  -  q'  sinh  36'  +  q'  sinh 


(3.13) 


(3.14) 


T2(8’)  =  1  -  2q *  cosh  2 S '  +  2q 


,4 


cosh  4b'  -  . .  . 


(3.15) 


T^U')  =  1  +  2q'  cosh  2 B '  ♦  2q 
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cosh  4b'  +  . . . 


(3.16) 
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(3.17) 


p  C 

T^(S')  =  cosh  8'  +  q'  cosh  36'  +  q'  cosh  58'  +  ... 


T02  "  T2(0)  =  1  "  2^'  +  ^  (3.18) 

TQ3  e  T3(0)  =  1  +  2q'  +  2q’4  +  ...  (3.19) 

Tq4  e  T4(0)  =  1  +  q'2  +  q'6  +  ...  (3.20) 

S  =  1  +  8q'2  -  8q'4  +  ...  (3.21) 


If  | 0 |  >  K  ,  the  elliptic  functions  can  be  calculated  from  the 
expressions  for  | 6 |  <  K  through  the  following  relations: 

cn(6  +  2nK;<)  =  (-1  )ncn(0;ic)  (3.22) 

sn(0  +■  2nK;<)  =  (-1)nsn(0;<)  (3.23) 

dn(0  +  2nK;tc)  =  dn(0;<)  (3.24) 

where  n  =  1,2,3... 
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The  development  and  presentation  of  the  expressions  for  the 
calculation  of  second-order  cnoidal  waves  has  been  completed  (Chapters 
II,  III,  and  Appendix  A).  The  next  chapter  will  discuss  the  shoaling 
and  refraction  of  second-order  cnoidal  waves  which  propagate  over  an 
irregular  bottom. 


CHAPTER  IV 


SHOALING  AND  REFRACTION  OVER  AN  IRREGULAR  SEA  BOTTOM 

Once  the  wave  height,  H,  wave  period,  T,  and  water  depth,  D,  are 
specified  at  a  given  location,  equations  presented  in  Chapters  II  and 
III,  and  in  Appendices  A  and  B  can  be  used  to  calculate  the  surface 
profile,  length,  celerity,  water  particle  velocities,  energy,  energy 
flux,  and  grou^  velocity  for  both  first-order  and  second-order  cnoidal 
waves  (provided  that  U  >  10).  The  direction  of  wave  propagation,  a  ,  is 
a  fourth  important  variable.  The  value  of  a  enters  into  the  calculation 
of  H,  as  well  as  other  important  engineering  quantities  such  as  wave- 
induced  currents  and  the  location  of  wave  breaking.  The  values  of  the 
four  variables,  H,  T,  D,  and  a  ,  are  necessary  at  every  location  where 
wave  and  wave- induced  current  calculations  are  needed. 

The  quantities  T  and  D  are  determined  or  specified  externally  to 
the  model.  Depths  are  specified  at  each  grid  node  prior  to  simulation. 
For  monochromatic  wave  modeling,  the  wave  period  is  conserved  as  the 
wave  transforms  over  the  irregular  bottom.  Therefore,  the  period 
specified  for  each  steady-state  simulation  is  the  period  used  for 
calculation  at  each  grid  node.  Consequently,  of  the  four  main 
variables,  only  H  and  a  remain  to  be  determined.  Their  calculation  is 
the  primary  function  of  the  numerical  model. 

Refraction  and  shoaling  (discussed  later  in  this  chapter)  are  the 
transformation  mechanisms  which  will  alter  H  and  a  as  the  wave 
propagates  over  the  irregular  bathymetry.  Energy  sources  and  sinks  are 
not  included  in  the  model.  The  action  of  wind  at  the  water-air  surface 
i3  therefore  ignored.  The  model  "assumes"  that  the  waves  are  formed 
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outside  the  modeled  area,  and  the  local  winds  do  not  add  appreciable 
energy.  The  effect  of  friction  at  the  bottom  is  also  omitted,  so  the 
model  will  give  best  results  if  used  for  coasts  that  have  sandy  rather 
than  muddy  bottoms.  (Muddy  bottoms  can  cause  significant  dampening 
effects.)  Dissipation  through  bottom  friction  could  be  straightforwardly 
included  in  the  model,  and  might  be  considered  in  future  developments. 

The  model  is  not  intended  for  use  in  the  surf  zone,  so  losses  from  wave 
breaking  are  not  considered.  Diffraction  caused  either  by  structures  or 
from  spatial  variation  of  bathymetry  is  not  included.  Therefore,  the 
model  is  an  "open  coast  model"  and  cannot  be  used  to  determine  wave 
heights  in  the  lee  of  structures,  or  in  locations  where  the  geometry  of 
the  bottom  appreciably  focuses  or  scatters  wave  energy.  The  addition  of 
any  such  transformations  will  demand  research  beyond  the  scope  of  this 
report. 

The  capability  to  calculate  heights  and  angles  of  cnoidal  waves 
directly  on  a  numerical  grid  is  an  important  feature  of  the  numerical 
model  developed  in  this  report.  No  other  known  finite-amplitude  wave 
transformation  model  has  this  capability.  All  other  known  finite- 
amplitude  models  and  many  of  the  small-amplitude  wave  transformation 
models  use  ray  theory,  a  classical  method  for  calculating  the  effects  of 
wave  refraction  over  an  irregular  bottom  (Munk  and  Authur,  1952). 

Defining  a  wave  front  as  a  line  of  constant  wave  phase  (such  as  the 
crest),  a  wave  ray  is  a  line  which  is  normal  to  the  wave  front.  In  the 
ray  tracing  method,  a  ray  equation  is  used  to  trace  rays  from  a  location 
where  wave  properties  are  known  across  a  region  where  wave  information 
is  desired.  A  second  equation  called  the  ray  separation  equation 
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determines  the  wave  height  resulting  from  refraction  by  assuming 
constant  energy  flux  between  adjacent  wave  rays. 

The  ray  tracing  technique  does  not  calculate  directly  on  a 
numerical  grid.  In  order  for  depth  information  to  be  available  at  the 
ray  location  and  for  the  calculated  heights  and  angles  to  be  given  at 
grid  nodes,  interpolation  is  necessary.  Therefore,  a  technique  first 
applied  by  Noda  et  al.  (1974)  was  adopted  in  this  study.  This  method 
calculates,  directly  on  a  numerical  grid,  the  wave  angles  and  heights 
resulting  from  the  shoaling  and  refraction  of  waves  over  an  irregular 
bottom.  Two  equations  are  used  to  predict  the  effects  of  shoaling  and 
refraction:  an  equation  for  the  calculation  of  the  angle  of  wave 
propagation,  and  a  conservation  of  energy  flux  equation,  which  is  used 
to  determine  wave  height. 

In  the  next  two  sections,  refraction  is  discussed,  and  the  wave 
angle  equation  is  derived.  Next,  the  mechanism  of  shoaling  is 
discussed,  and  the  energy  conservation  equation  used  to  predict  wave 
height  will  be  derived.  Finally,  expressions  for  energy  flux  and 
average  wave  energy  per  unit  surface  area  of  second-order  cnoidal  waves 
will  be  presented. 

4.1  Wave  Refraction 

As  can  be  seen  from  Equation  2.65,  at  a  first  approximation,  the 
wave  celerity  varies  as  the  square  root  of  the  depth.  Therefore,  if  the 
depth  is  not  constant  along  the  wave  front,  the  local  wave  celerity  will 
also  vary  along  the  wave  front.  With  different  portions  of  the  wave 
front  traveling  at  different  speeds,  the  wave  begins  to  align  parallel 


to  the  bottom  contours.  This  transformation,  called  refraction  in  a 
direct  analogy  with  optics,  can  significantly  affect  the  wave  height  by 
causing  the  convergence  or  divergence  of  wave  energy. 

Refraction  is  especially  important  for  the  calculation  of  sediment 
transport.  The  magnitude  of  wave-induced  longshore  current,  the 
principal  cause  of  longshore  sediment  transport  and  beach  change,  is  a 
function  of  the  angle  of  wave  propagation  at  wave  breaking  (SPM,  1984, 
p.  4-54).  Also,  since  the  position  of  wave  breaking  is  a  function  of 
wave  height  and  water  depth,  the  location  of  the  wave  breaking  point  is 
controlled  in  part  by  wave  refraction.  The  location  of  wave  breaking 
determines  the  width  of  the  area  over  which  the  longshore  current  and 
intense  sediment  transport  occur. 

Traditionally,  coastal  engineers  and  scientists  have  broadened  the 
usage  of  the  terminology  "wave  refraction"  to  include  wave  shoaling  in 
addition  to  pure  refraction.  Thus,  a  numerical  model  of  wave  refraction 
usually  implies  that  both  shoaling  and  refraction  have  been  considered. 
All  three  of  the  finite-amplitude  refraction  models  for  an  irregular 
bathymetry  discussed  in  the  introduction  (Chu,  1975;  Headland  and  Chu, 
1984;  and  Oh  and  Grosch,  1985),  predicted  that  finite-amplitude  waves 
refract  less  than  do  small-amplitude  waves.  Skovgaard  and  Petersen 
(1977)  calculated  the  refraction  of  first-order  cnoidal  waves  for 
straight  and  parallel  bottom  contours  with  the  same  result. 

Wang  and  Le  Me'haute  (1980)  found  that  the  prediction  of  less 
refraction  for  finite-amplitude  waves  did  not  match  well  to  one  set  of 
laboratory  experiment  observations.  They  developed  a  "hybrid"  model, 
which  uses  the  assumption  of  a  plane  bottom.  The  hybrid  model  employs 


onoidal  wave  theory  for  shoaling  calculations  and  small -amplitude  theory 
for  refraction  calculations.  This  combination  matched  the  experimental 
data  for  waves  of  relatively  high  deep-water  steepness.  The  results 
showed  that  both  the  data  and  the  hybrid  model  indicated  less  refraction 
than  small -amplitude  theory,  but  more  refraction  than  predicted  by 
finite-amplitude  theory.  However,  the  differences  in  wave  refraction 
which  motivated  the  hybrid  model  were  observed  in  the  laboratory 
experiment  and  calculated  in  the  numerical  model  at  wave  breaking,  where 
the  numerical  model  must  estimate  the  breaking  point,  as  well  as 
calculate  refraction.  Therefore,  it  is  difficult  to  determine  if  the 
differences  between  calculated  and  observed  wave  angle  were  due  to 
errors  in  the  physical  model  observations,  under-prediction  in  the 
numerical  modeling  of  refraction  by  finite-amplitude  theory,  difference 
between  the  breaking  criterion  used  in  the  numerical  model  and  the 
breaking  position  observed  in  the  physical  modeling,  or  some  combination 
of  the  three.  The  question  of  whether  finite-amplitude  theory  under¬ 
predicts  refraction  cannot  be  answered  without  more  experimental 
results.  The  hybrid  model  seems  to  be  too  drastic  a  step  based  on  the 
observations  of  one  experiment.  Additional  well-controlled  laboratory 
experiments  are  required  to  resolve  this  point. 

4.2  Derivation  of  the  Wave  Angle  Equation 

The  equation  and  procedure  for  calculating  the  wave  angle  as 
introduced  by  Noda  et  al.  (1974)  was  derived  for  small-amplitude  wave 
theory.  It  is  based  on  the  irrotat tonality  of  wave  number.  The  wave 
number  in  small-amplitude  theory  is  defined  as  the  gradient  of  the  phase 


function  of  the  surface  profile.  Since  the  phase  functions  of  cnoidal 
and  small -amplitude  waves  are  different,  the  wave  angle  equation  must  be 
rederived  for  use  with  cnoidal  wave  theory.  The  following  development 
of  the  wave  angle  equation  is  patterned  after  a  derivation  for  small- 
amplitude  theory  found  in  Dean  and  Dalrymple  (1984,  Chapter  4). 

The  surface  profile  of  a  second-order  cnoidal  wave  (Equation  2.66) 
is  repeated  in  Equation  4.1  for  easy  reference. 

n  =  D  +  A1  cn2  0  +  A2  on**  0^  (4.1) 

For  a  wave  which  is  traveling  in  an  arbitary  direction  in  the  x-y  plane, 

9  is  given  as 

(4.2) 

where  L/cos  a  and  L/sin  a  ,  the  projections  of  the  wave  crest  on  the  x 

and  y  axis,  respectively;  and  a  ,  the  direction  of  wave  propagation,  are 

shown  in  F iqure  4.1. 

The  phase  function  can  be  written  as 

9  =  2K a  (4.3) 


0  =  2K  I  ^  cos  a  +  ^  sin  a  -  ^ 


I 

Fiqure  4.1.  Definition  sketch  for  a  ,  L/cos  a  ,  and  L/sin  a  . 

2  4 

As  illustrated  in  Appendix  A,  the  periodic  functions  cn  e  and  cn  e 
of  Equation  4.1  have  a  non-constant  period,  2K.  However,  the  crest  of 
the  wave  of  Equation  4.1  occurs  for  ft  =  n;  n  =  0,1,2,....  ,  regardless 
of  the  value  of  K.  Therefore,  the  normal  vector  to  the  crest  (and  thus 
the  direction  of  wave  propagation)  can  be  determined  by  taking  the 
gradient  of  Q  ,  as  in 

i  +  (^t^J  J  (4-5) 

where  1  and  ]  are  the  unit  vectors  in  the  x  and  y  directions. 

From  a  vector  identity  that  states  that  the  curl  of  the  gradient 
of  a  scaler  is  equal  to  zero, 


VQ 


(cos  q 
L  , 
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V  x  VQ  =  0 


(4.6) 


which  can  be  rewritten  as 


This  is  the  governing  equation  used  in  the  numerical  model  to 
calculate  a  on  a  fixed  grid.  The  finite  difference  representation  of 
Equation  4.7  will  be  presented  in  Chapter  V.  It  is  important  to  note 
that  wavelength,  L,  is  a  function  of  wave  height,  H,  in  a  finite- 
amplitude  wave  theory.  This  can  be  shown  using  the  basic  relationship, 
L  =  CT,  and  from  Equation  2.65,  where  it  is  seen  that  wave  celerity,  C, 
(from  e  =  H/D)  is  a  function  of  H.  This  has  important  implications  for 
the  calculation  of  refraction  based  on  a  finite-amplitude  wave  theory. 
First,  refraction  causes  changes  in  wave  height,  which  will  in  turn 
affect  the  calculation  of  wave  angle  in  Equation  4.7.  Shoaling  also 
causes  changes  in  wave  height,  as  will  be  discussed  below.  Therefore, 
the  dependence  of  refraction  on  wave  height  couples  shoaling  and 
refraction. 

The  equation  derived  for  the  calculation  of  the  angle  of  wave 
propagation  (Equation  4.7)  reduces  to  a  simpler  form  often  used  for 
refraction  studies.  Snell's  law  (Equation  4.8)  is  used  to  study 
refraction  in  analytical  and  theoretical  studies.  It  is  also  commonly 
used  if  the  assumption  of  straight  and  parallel  bottom  contours  is 
Justified  or  used  locally  over  short  distances  for  a  more  complicated 
bathymetry. 


sin  a  .  . 

— ^ —  =  constant 


(4.8) 
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Establishing  a  coordinate  system  so  that  the  x-axis  is  normal  to  shore 
and  the  y-axis  is  along  the  shore,  the  y-derivatives  vanish  for  the  case 
of  straight  and  parallel  bottom  contours.  Equation  4.7  becomes 


This  equation  can  be  converted  to  Snell's  law  (Equation  4.8)  by 
multiplying  by  T,  using  the  relationship,  C  =  L/T,  and  integrating  in 
the  x  direction.  The  constant  in  Equation  4.8  is  usually  evaluated  as 
sin  aQ/Lo  in  which  the  subscript  "o"  denotes  known  quantities  at  a 
reference  depth,  usually  in  deep  water. 

4.3  Wave  Shoaling 

Changes  in  water  depth  can  also  cause  changes  in  wave  height  due 
to  the  mechanism  known  as  wave  shoaling.  As  a  wave  "shoals,"  the  wave 
height  adjusts  to  compensate  for  variations  in  the  velocity  of  wave 
energy  propagation  caused  by  changes  in  depth.  In  general,  as  a  wave 
first  begins  to  shoal  there  is  a  region  where  the  wave  height  will 
decrease  slightly.  Then,  as  the  wave  enters  progressively  shallower 
water,  the  wave  height  will  increase.  Observations  of  waves  just  before 
the  breaking  point  show  that  a  rapid  increase  in  height  can  occur.  The 
investigation  of  the  equations  governing  energy  propagation  of  second- 
order  cnoidal  waves  will  exhibit  the  effect  of  depth  on  shoaling. 

The  average  rate  of  wave  energy  propagation  through  a  vertical 
cross  section  normal  to  the  direction  of  wave  propagation,  called  the 
energy  flux,  F,  is  defined  in  Equation  4.10  (Phillips,  1977,  p.  25).  (An 


equation  for  the  calculation  of  energy  flux  for  second-order  cnoidal 
waves  is  derived  in  Appendix  B  and  is  presented  and  discussed  at  the  end 
of  this  chapter.)  The  quantity  E  is  the  average  energy  per  unit  surface 


area  of  the  wave,  and  the  group  velocity,  Cg,  is  the  velocity  of  wave 
energy  propagation. 


f  =  E  2  (4.10) 

g 

Both  group  velocity  and  energy  flux  are  vector  quantities,  but  the 
vector  notation  will  be  omitted  for  convenience,  except  if  the 
directionality  of  the  two  quantities  is  emphasized.  The  customary  term, 
group  "velocity,"  will  be  used  instead  of  the  term,  group  "speed,"  when 
referring  to  C  .  Equations  4.11  and  4.12,  for  C  and  E,  are  derived  in 

o  © 

Appendix  B  (Equations  B.88  and  B.83).  EQ,  E1*  and  F1*  functions  of 
\  and  u  ,  are  defined  in  Equations  B.83  and  B.46. 


C 

g 


E  =  pgH2  (EQ  +  eE}) 


(4.11) 


(4.12) 


At  a  first  approximation,  Equation  4.11  indicates  that  C  varies  with 

O 

the  square  root  of  depth.  As  will  be  discussed  below,  it  is  assumed 

that  F  is  conserved  as  the  wave  shoals.  Therefore,  from  Equation  4.10, 

if  C  decreases  with  decreasing  depth,  then  E,  in  order  to  conserve  F, 
© 

p 

must  compensate  by  increasing.  Since  E  is  proportional  to  Hc,  changes 
in  E  cause  corresponding  changes  in  H. 


There  have  been  numerous  investigations  devoted  to  the  subject  of 
shoaling  of  finite-amplitude  waves.  All  the  researchers  have  found  that 
predictions  of  wave  shoaling  using  small -amplitude  theory  underestimate 
shoaling  rates  predicted  by  either  finite-amplitude  theories  or 
exhibited  in  experimental  data.  Le  Me'haute  and  Webb  (1964)  used  the 
conservation  of  energy  flux  to  study  the  shoaling  of  third-order  Stokes 
waves.  Koh  and  Le  Mehaute  (1966)  found  that  fifth-order  Stokes  waves 
have  a  slightly  lower  shoaling  coefficient  than  do  third-order  Stokes 
waves.  An  equation  to  calculate  shoaling  for  first-order  cnoidal  waves 
was  derived  by  Svendsen  and  Brink-KJaer  (1972).  They  supplied  tables 
which  can  be  used  to  predict  wave  height  and  wave  length  in  shallow 
water  using  deepwater  inputs.  Shuto  (1974)  derived  a  law  of  shoaling 
from  basic  hydrodynamic  equations,  and  presented  expressions  for  the 
engineering  calculation  of  shoaling  in  three  different  regions  defined 
by  the  Ursell  number.  The  expression  given  for  the  region  U  >  30 
corresponds  to  first-order  cnoidal  theory.  Sakai  and  Battjes  (1980) 
calculated  shoaling  based  on  Cokelet's  high-order  and  essentially 
mathematically  exact  wave  theory,  and  compared  the  results  with  those  of 
the  above  researchers.  They  found  that  Cokelet's  numerical  theory 
predicted  slightly  higher  shoaling  than  the  above-mentioned  studies. 
Walker  and  Headland  (1982)  compared  several  of  the  above  theoretical 
shoaling  calculations  with  empirically  based  wave  breaking  indices.  They 
found  that  the  theoretical  results  predicted  greater  shoaling  rates  than 
indicated  by  the  empirical  breaking  indices.  However,  they  concluded 
that  the  wave  flume  data  on  which  the  empirical  curves  are  based  may  be 
contaminated  by  unwanted  harmonics  resulting  from  sinusoidal  motion  of 


the  wave  generator.  Isobe  (1985)  used  small-amplitude  theory  (for 
deeper  water)  and  first-order  cnoidal  theory  (for  shallower  water)  to 
compute  shoaling,  and  compared  results  to  experimental  data.  The 
comparison  was  favorable  for  waves  of  small  to  intermediate  deepwater 
wave  steepness. 

4.4  Derivation  of  the  Conservation  of  Energy  Flux  Equation 

The  conservation  of  energy  principle  applied  in  a  control  volume 
approach  shows  that 

V  •  P  =  0  (4.13) 

under  the  following  assumptions: 

a.  The  bottom  slope  is  sufficiently  mild  so  that  reflection  is 
negligible. 

b.  There  are  no  energy  sources  or  sinks. 

c.  The  solution  is  for  the  steady  state. 

Equation  4.13  is  the  basis  for  the  calculation  of  wave  height  in  the 
numerical  model.  The  technique  and  the  finite  difference  equations  will 
be  presented  in  Chapter  V.  Use  of  this  equation  in  the  numerical  model 
requires  an  expression  for  the  energy  flux  of  second-order  cnoidal 
waves.  The  development  of  this  expression  was  a  nontrivial  task  and  is 
a  major  accomplishment  of  this  study.  The  derivation  of  an  expression 
for  the  energy  flux  of  second-order  cnoidal  waves  is  included  in 
Appendix  B,  and  the  resulting  equation  for  F  (Equation  B.46),  is 
repeated  as 


F  =  pgH2  yj  gD  (  Fq  +  eF ^ ) 


=  pgH2  yj  gD  |^(-X  +  2y  +  4Xy  -  X2  -  3u2)  ^ 
+  e  (-4X  +  8 y  +  53Xy  -  1 2 X2  -  60y2 


+  53X2y  -  120Xy2 


8x: 
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(4.14) 


For  comparison,  only  one  other  source  could  be  found  which  gave 
expressions  for  F  and  E  for  second-order  cnoidal  waves.  The  equation 
for  F  (Sarpkaya  and  Isaacson,  1981,  Chapter  4),  derived  from  the  second- 
order  cnoidal  theory  presented  by  Laitone  (I960)  is 


Fs&i  =  pgI)2  V^"e‘ 
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(4.15) 
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From  X  =  k'/k  (Equation  2.70),  and  from  y  =  Eg/K  K  (Equation  2.71) 
Equation  4.14  can  be  shown  to  be  equivient  to  Equation  4.15.  This 
result  is  expected  since  both  Laitone's  derivation  of  second-order 
cnoidal  wave  theory  and  the  derivation  used  in  this  thesis  (Isobe  and 
Kraus,  1983b)  use  Stokes'  second  definition  of  wave  celerity  (Equation 
2.21).  The  expression  for  E  (Equation  B.83)  was  also  found  to  be 
equivalent  to  that  given  by  Sarpkaya  and  Isaacson  ( 1 98 1 ,  Chapter  4).  A 


i  similar  check  for  C  Equation  B.88  is  not  possible  since  no  other 

t  © 
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expression  for  the  group  velocity  of  second-order  cnoidal  waves  is 
known . 

The  basis  has  been  established  for  the  calculation  of  the  shoaling 
and  refraction  of  second-order  cnoidal  waves  over  an  irregular 
bathymetry.  The  next  chapter  describes  the  numerical  model  designed  to 
execute  this  task. 


CHAPTER  V 


DESCRIPTION  OF  THE  NUMERICAL  MODEL 
5.1  Overview  of  Technique 

i 

j  As  was  discussed  in  Chapter  IV,  shoaling  and  refraction  result 

from  changes  in  group  velocity,  Cg  and  wave  celerity,  C,  respectively, 
as  waves  propagate  over  varying  water  depth.  For  small-amplitude  wave 
theory,  C  and  C  are  functions  of  wave  period  and  water  depth  only. 

O 

Since  neither  the  wave  period  nor  water  depth  are  affected  by  shoaling 
or  refraction,  shoaling  and  refraction  can  be  directly  and  independently 
calculated  in  a  model  based  on  small-amplitude  wave  theory.  For 
instance,  the  wave  height  at  any  location  can  be  calculated  by 
determining  a  shoaling  coefficient,  Kg,  and  a  refraction  coefficient, 

J  Kr,  each  calculated  independently  of  the  other.  With  the  deep  water 

f 

wave  height  signified  by  HQ,  the  wave  height  from  a  small -amplitude 
model  is  given  as  H  =  Kr^s^o' 

|  For  finite-amplitude  wave  theory  (in  this  report,  both  first  and 

second-order  cnoidal  theory),  both  Cg  (Equation  B.88)  and  C  (Equation 
2.65),  are  functions  of  wave  height,  as  well  as  wave  period  and  water 
depth.  Therefore,  since  both  shoaling  and  refraction  cause  changes  in 
wave  height,  there  is  an  interdependence  between  the  two  transfor¬ 
mations.  Consequently,  the  calculation  of  wave  height  caused  by 
shoaling  and  refraction  using  a  finite-amplitude  theory  is  not  as  direct 
as  with  small-ampl itude  theory.  In  the  numerical  model  developed  in 
this  study,  an  iterative  technique  is  used  to  account  for  the  inter- 

i 

;  dependence  of  shoaling  and  refraction.  First  an  estimate  for  the  angle 

t 

of  wave  propagation  is  calculated,  followed  by  the  calculation  of  an 
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estimate  for  the  wave  height.  The  process  is  repeated  until  the  change 
in  wave  height  from  one  iteration  to  the  next  is  below  a  specified 
tolerance.  The  number  of  iterations  varies  from  two  to  six  for  con¬ 
vergence  to  a  tolerance  that  is  about  0.1%  of  the  wave  height.  In 
general,  more  iterations  are  necessary  near  to  the  wave  breaking  point, 
where  the  wave  height  is  changing  rapidly,  than  are  needed  near  to  the 
seaward  boundary  of  the  model. 

Figure  5.1  shows  the  arrangement  of  the  numerical  grid.  The  grid 
is  aligned  so  that  the  y-axis  (J=1,n)  lies  in  the  alongshore  direction 
and  the  x-axis  (i=1,m)  lies  normal  to  the  shoreline.  With  the  origin  of 
the  grid  located  on  dry  land  it  is  easier  to  adjust  the  position  of  the 
seaward  boundary  of  the  area  to  be  modeled.  The  x-axis  is  placed  so  that 
all  the  area  to  be  modeled  has  positive  y-coordinates.  For  ease  of 


description  in  the  remainder  of  this  chapter,  the  term  "row"  is  used  to 
mean  a  line  of  constant  x  (roughly  parallel  to  the  shoreline),  where  the 
nodes  of  the  numerical  grid  are  designated  (j,i);  J=1,n  ;  i=constant. 
(The  order  of  all  field  arrays  in  the  computer  code  is  (J,i)  to  allow 
for  more  efficient  use  of  computer  memory.) 

Figure  5.1  also  contains  a  definition  sketch  for  the  wave  angle. 
Both  input  and  output  wave  angles  are  measured  counterclockwise 
(direction  from  which  the  waves  are  incident)  from  a  line  parallel  to 
the  x-axis.  For  internal  model  calculations  180°  is  added  to  these 
defined  angles  so  that  calculations  can  use  a  definition  of  angle  which 
matches  the  sign  of  the  trigonometric  functions  and  the  coordinate 
system  of  the  numerical  grid.  For  the  internal  calculations  wave  angle 
is  measured  counterclockwise  (direction  which  the  wave  is  traveling) 
from  a  line  parallel  to  the  x-axis. 

The  wave  parameters  of  height  (H),  period  (T),  and  angle  (a)  are 
the  basic  inputs  to  the  model.  The  numerical  model  calculates  for  a 
steady  state  condition,  so  each  wave  condition  requires  a  separate 
simulation.  Each  simulation  assumes  that  wave  conditions  are  constant 
for  the  time  that  it  would  take  a  wave  to  travel  from  the  seaward 
boundary  to  the  breaking  point.  The  wave  period  is  specified  as  a 
constant  for  all  grid  nodes  during  each  simulation.  Wave  height  and 
wave  angle  can  be  specified  as  constant  along  the  whole  seaward  boundary 
or  they  can  vary  for  each  grid  node  on  the  seaward  boundary,  resulting 
from  the  calculations  of  a  companion  model  used  for  deeper  water. 
(Methods  of  connection  with  a  companion  model  will  not  be  investigated 
in  this  thesis.)  The  boundary  conditions  on  the  lateral  sides  of  the 


grid  are  the  simplest  possible.  The  lateral  boundary  values  of  H,  a, 
and  the  energy  flux,  F,  take  the  values  of  the  node  next  to  the 
boundary.  Tests  showed  that  the  model  is  insensitive  to  the  lateral 
boundary  conditions.  As  long  as  the  lateral  boundaries  are  not  close  to 
areas  of  interest,  the  model  results  are  negligibly  affected  by  the 
lateral  boundary  conditions.  In  the  present  version  of  the  model  the 
lateral  boundaries  are  assumed  to  be  in  water;  therefore,  the  model  is 
configured  for  simulations  on  a  relatively  straight  coastline  with  no 
significant  embayments. 

Since  the  computational  scheme  marches  toward  the  shoreline,  a 
boundary  condition  along  the  shoreline  is  not  essential  from  a  numerical 
standpoint.  However,  since  the  model  is  not  intended  for  use  shoreward 
of  the  wave  breaking  point,  the  computations  are  halted  before  the 
marching  scheme  reaches  the  shoreline.  For  operation  in  a  plane  beach 
mode,  the  wave  heights  will  be  identical  along  each  row.  In  thi3  case  a 
standard  breaking  criterion,  Hb  =  0.8D,  is  used.  The  simulation  is 
halted  if  the  breaking  criterion  is  exceeded.  For  the  case  of  irregular 
bathymetry,  wave  heights  will  vary  along  each  row,  and  the  onset  wave 
breaking  will  not  occur  at  the  same  row  throughout  the  modeled  area.  If 
the  calculated  wave  height  exceeds  the  breaking  criterion  at  a  node,  the 
wave  height  is  truncated  so  that  H  =  0.8D  .  This  technique  allows  the 
location  of  wave  breaking  to  be  determined  for  the  entire  modeled  area. 
The  simulation  halts  if  the  breaking  criterion  has  been  exceeded  at  all 
nodes  on  a  row. 

Judgement  must  be  exercised  in  applying  model  results  near  the 
breaking  point.  The  perturbation  solution  for  the  second-order  cnoidal 
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wave  theory  used  in  this  report  is  based  on  the  assumption  that  the 
i  ratio  of  wave  height  to  water  depth  is  small;  therefore,  model  results 

as  H/D  approaches  1.0  are  suspect.  Futhermore,  as  a  wave  nears 
breaking,  the  assumption  of  permanent  form,  which  is  fundamental  to  the 
derivation  of  cnoidal  wave  theory,  is  violated.  This  difficulty  is  not 
I  unique  to  this  numerical  model;  it  occurs  in  all  shallow  water  wave 

transformation  models  currently  available  for  engineering  use. 

A  flow  chart  of  the  numerical  model  is  presented  in  Figure  5.2. 

One  of  chief  assets  of  this  model  is  its  simplicity.  The  difficult  work 
has  been  the  development  of  the  equations  for  the  calculation  of  shoaling 
and  refraction  directly  at  grid  points  (Chapter  IV  and  Appendix  B),  and 
the  development  of  expressions  for  the  efficient  calculation  of  elliptic 
integrals  and  Jacobian  elliptic  functions  (Appendix  A).  This  is  the 
advantage  of  using  an  analytical  wave  theory  instead  of  a  numerical  wave 
theory.  The  model,  written  in  FORTRAN  V,  is  divided  into  three  main 
segments,  an  input  segment,  the  main  calculation  segment,  and  an  output 
segment. 

5.2  Model  Output  and  Input 

The  primary  results  of  the  model  are  wave  height  and  wave  angle  at 
each  node  on  the  numerical  grid.  Depth,  wave  length,  wave  celerity, 
group  velocity,  and  Ursell  number  can  also  be  included  in  the  output  by 
specifying  a  keyed  number  in  the  model  input  file.  The  output  file 
consists  of  one  grid  for  each  of  the  chosen  quantities  with  the  values 
of  the  particular  quantity  given  at  each  node. 

The  input  file  read  into  the  model  at  the  beginning  of  each 
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Figure  5.2.  Flow  chart  of  numerical  model. 
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simulation  contains  values  of  the  acceleration  of  gravity,  the  density 
of  water,  the  desired  convergence  criterion,  a  correction  to  the  water 
depths  to  account  for  changes  in  water  depth  due  to  tide  or  storm  surge, 
the  output  specifier,  a  flag  which  triggers  the  input  or  calculation  of 
the  depth  grid,  and  constants  for  use  in  calculating  an  idealized 
bathymetry.  Also  read  in  from  the  input  file  are  the  grid  spacing  in 
the  both  the  x  and  y  directions.  The  numerical  grid  is  a  fixed 
rectangular  grid,  and  the  values  of  the  grid  spacing,  Ax  and  Ay,  do  not 
have  to  be  equal.  Also  read  from  the  input  file  is  a  flag,  through 
which  the  user  can  specify  that  either  first  or  second-order  cnoidal 
wave  theory  be  used. 

The  depths  at  each  grid  node  are  read  from  the  input  file  for  the 
case  of  prototype  bathymetry  or  can  be  calculated  for  the  case  of 
idealized  bathymetry  using 


D  .  =  a  +  bx  -  (d)(i/m)  cos  [2n(J/n)]  +  (d)(i/m)  (5.1) 

J  > 1 

This  equation  is  composed  of  a  plane  beach  component  and  a  longshore 
rythmic  component.  Such  a  topography  is  an  idealization  of  periodic  rip 
channels  superimposed  on  a  uniformly  sloping  bottom,  as  is  frequently 
seen  along  sandy  coasts.  The  symbols  for  the  plane  bottom  component, 

"a"  and  b,  are  the  depth  at  i=0  and  the  bottom  slope,  respectively.  The 
rythmic  component  varies  the  depth  in  the  longshore  direction  with  a 
sinosoidal  pattern  which  has  a  phase  determined  by  the  normalized 
alongshore  position.  The  sinisoidal  pattern  has  an  maximum  amplitude  of 
2d  at  the  seaward  boundary  and  decays  towards  the  shoreline  depending 
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upon  the  value  of  the  normalized  onshore-offshore  position.  If  d  is 
positive,  a  trench  is  created.  If  d  is  negative,  a  shoal  is  created. 

5.3  Calculating  Wave  Height  and  Wave  Angle 

Four  subroutines  are  used  for  the  main  computational  scheme, 

ELLIP,  LENGTH,  ANGLE,  and  EFLUX.  Before  calculations  can  begin  for  the 

interior  of  the  grid,  wave  celerity  (C),  wave  length  (L),  and  energy 

flux  (F)  must  be  calculated  on  the  seaward  boundary.  Since  C,  L,  and  F 

are  all  functions  of  X  and  u  (Equations  2.65,  L  =  CT,  and  4.14),  ELLIP 

is  used  to  calculate  X  and  y.  Then  C  and  L  are  calculated  in  LENGTH. 

Finally  the  calculation  of  F  in  EFLUX  completes  the  specification  of  the 

seaward  boundary  conditions.  (The  wave  angle,  a  is  read  from  the  input 

file.)  Once  the  above  quantities  are  calculated  on  the  seaward  boundary 

the  computational  scheme  marches  toward  shore  (in  the  negative  x 

direction)  one  row  at  a  time.  The  numerical  scheme  is  explicit,  with 

information  from  from  ith  row  used  to  calculate  values  at  the  (i-1)th 

row.  Calculations  are  repeated  until  H  converges  to  within  the 

specified  tolerance  at  each  node  on  the  row  before  preceeding  to  the 

next  row.  Before  the  first  iteration  on  each  row,  as  an  initial  guess, 

the  wave  height  at  each  node  is  set  equal  to  the  wave  height  of  the  node 

0  k 

from  the  previous  row  in  the  same  column,  (H  ) j  j  1  =  (H  )j  ^  ;  J  =  1 ,n, 
where  the  superscript,  k,  denotes  iteration  level.  The  four  subroutines 
used  in  the  calculations  are  described  in  the  following  sections. 

5.3.1  Subroutine  ELI. IP 

Since  the  first  and  second-order  cnoidal  wave  quantities  have  been 
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expressed  in  terms  of  X  and  y,  evaluation  of  X  and  y  is  necessary.  An 
equation  for  the  evaluation  of  X  derived  from  Equations  2.70,  3.6,  and 
3.7  is 


x 


(5.2) 


An  equation  for  the  evaluation  of  y  derived  from  Equations  2.71,  3.6, 
3.8,  and  3.9  is 


y  = 


(T02) 


_ ? _  _  S  +  (T 

(-lnq1 )  5  U03} 


(5.3) 


The  complementary  nome,  q'  is  the  basic  quantity  used  in  the  calculation 
of  these  two  quantities,  as  well  as  all  other  cnoidal  wave  quantities. 

An  equation  for  the  calculation  of  q'  can  be  derived  from  Equation  2.85, 
the  dispersion  relation  for  second-order  cnoidal  waves  repeated  here  as 


16  2„2 
— 2  K  K  = 


$HT_ 

D2 


1  +  e 


(M 


(5.4) 


If  the  second  term  in  the  brackets  of  the  above  equation  is  neglected, 
the  dispersion  relation  for  first-order  cnoidal  waves  results.  From  k 
and  K  (Equations  3.5  and  3.7),  Equation  5.4  can  be  rewritten  as 


4 

3 


(T02> 


(-lnq')2 


(5.5) 


Solving  this  equation  for  q'  gives 
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q'  = 

For  a  simulation  employing  first-order  cnoidal  wave  theory,  the  higher 
order  contribution  to  q'  (the  coefficient  of  e  in  Equation  5.6)  is  set 
equal  to  zero. 

Note  that  both  T Q2  (Equation  3.18)  and  X  (Equation  2.70)  are 
functions  of  q'.  Therefore,  an  iterative  solution  is  used  to  solve  for 
q'.  Since  the  iteration  for  q'  is  inside  the  iteration  loop  for  H,  a 
technique  using  under-relaxation  was  developed  to  converge  on  q'  in  as 
few  iterations  as  possible.  Without  the  relaxation  technique  it  can 
take  up  to  13  iterations  to  converge  to  (q,k  -  q,lc_1  )/q'k_1  <  0.001  for 
U  =  10.  This  would  increase  computer  time  since  q'  must  be  calculated 
at  each  node.  The  relaxation  equation  is 

q,k  =  q'*  -  (q,#  -  q'k_1>  R  (5.7) 

where  R  is  the  relaxation  parameter  and  q'  is  the  value  of  q'  at  the 
kth  iteration  before  relaxation.  Initial  tests  showed  that  a  constant 
value  for  R  does  not  produce  acceptable  results  for  the  full  range  of 
Ursell  numbers  applicable  to  cnoidal  wave  theory  (U  >  10).  Difficulties 
with  convergence  occur  at  smaller  Ursell  numbers  so  the  optimum  value  of 
R  was  determined  for  a  series  of  Ursell  numbers  between  10  and  100. 

Since  the  optimum  value  varies  not  only  with  the  Ursell  number  but  also 
with  the  value  of  the  perturbation  parameter,  representative  values  of  e 
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were  used,  and  an  exponential  curve  was  fit  through  the  optimum  R  values 
for  U  =  10  (e  =  0.1)  and  U  =  30  (e  =  0.2),  giving 

R(U>  =  e«-°-61-0.054U>  (5.8> 

The  calculation  for  q'  converges  to  (q'^  -  q'^"^)/q'^_^  <  0.001  after 
only  three  iterations  for  U  >  10,  using  Equation  5.8  to  determine  the 
relaxation  parameter.  For  the  interior  of  the  grid,  since  calculations 
are  repeated  at  least  twice  because  of  the  iteration  for  wave  height, 
only  two  iterations  for  q'  are  needed  during  each  iteration  for  wave 
height.  For  calculations  on  the  seaward  boundary  six  iterations  are 
used  to  ensure  very  accurate  determinations  of  q'  for  the  boundary 
conditions.  (This  row  is  not  included  inside  the  iteration  for  H.)  The 
model  determines  the  least  accurate  convergence  for  each  row  and  those 
values  are  printed  at  the  end  of  the  simulation.  The  method  of  using  a 
fixed  number  of  iterations  has  the  advantage  of  not  having  to  check  a 
convergence  criterion  with  FORTRAN  "IF"  statements. 

5.3.2.  Subroutine  LENGTH 

Once  X  and  u  are  calculated,  the  wave  celerity  is  calculated  from 
Equation  2.65.  The  zeroth-order  contribution  to  C,  i.e.,  CQ,  is  also 
calculated  for  use  in  the  calculation  of  Cg  and  F  in  Subroutine  EFLUX. 
The  wave  length,  L,  is  then  calculated  using  the  relation  L  =  CT. 

For  a  simulation  employing  first-order  cnoidal  wave  theory,  the 
second-order  contribution  to  C,  i.e.,  C2,  is  set  equal  to  zero. 
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5.3.3  Subroutine  ANGLE 


Equation  4.7,  the  governing  equation  for  the  calculation  of  the 
angle  of  wave  propagation  is  repeated  here  as 


3  /sin  q\  3  /cos  a\  _  Q 
3x\L  /  3y  \  L  / 


(5.9) 


The  FTCS  (forward  time,  central  space,  where  the  x  is  considered  "time") 
explicit  finite  difference  scheme  used  in  the  model  is  given  as 
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(5.10) 


This  equation  is  solved  for  a  .  ;  j  =  1,n,  (The  scheme  marches  in  the 

J  * 1_  • 

minus  x-direction . )  which  give; 


(5.11) 


In  their  numerical  model  of  small -amplitude  wave  refraction,  Perl in  and 
Dean  (1983)  found  that  an  implicit  version  of  this  finite  difference 
scheme  can  cause  oscillations  with  a  wavelength  equal  to  2Ay.  They  added 
a  "dissipative  interface"  to  the  finite  difference  scheme  to  control 
these  spurious  oscillations.  This  problem  has  not  yet  been  encountered 
in  the  model  developed  in  this  report.  Equation  5.11  has  been  well- 
behaved,  possibly,  because  only  smoothly  varying  idealized  bathymetries 


have  been  used  in  the  simulations.  Future  testing  of  the  model  will 
involve  optimizing  the  finite  difference  technique. 

5.3.4  Subroutine  EFFLUX 

The  governing  equation  for  the  calculation  of  energy  flux 
(Equation  4.13)  is  repeated  here  as 

7  •  ~F  =  0  (5.12) 


Expanding  this  equation  gives 


— ^  (F  cos  a)  +  — ^ -  (F  sin  a)  =0 
3x  3y 


(5.13) 


Again  as  in  the  finite  differencing  of  the  wave  equation  an  explicit 
FTCS  finite  difference  algorithm  is  used.  This  results  in 


F,  , 
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J , i-1  (cos  “)j>i_1 


(Fcos  «)Jit  ♦  2jf[(Fsin  a)Jitji  -  (Fsin 


From  Equation  4.14,  F  can  be  rewritten  as 


(5.14) 


fj,i  =  hj,i  F'.i 


(5.15) 


where  F  is  given  as 
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The  wave  height  is  calculated  using 


F  .  \1/2 

-r1)  (5.i7 

FJ,i/ 

where  FQ  and  F1  are  defined  in  Equation  B.H6. 

Finally,  C  is  calculated  from  Equation  B.88,  after  a  value  for  E( 

O 

is  determined  using  Equation  B.83. 

If  first-order  cnoidal  wave  theory  is  selected  in  the  input  file, 
the  values  of  F1  and  E1  are  set  equal  to  zero. 

The  next  chapter  will  present  results  from  simulations  of  the 
shoaling  and  refraction  of  both  first  and  second-order  cnoidal  waves. 


CHAPTER  VI 

NUMERICAL  MODEL  RESULTS 

Results  from  numerical  simulations  of  the  shoaling  and  refraction 
of  shallow  water  waves  will  be  presented  in  this  chapter.  Section  A.1 
and  Section  A. 2  contain  results  for  shoaling  and  refraction,  respec¬ 
tively,  over  a  plane  bottom.  In  these  sections,  a  comparison  is  made 
among  the  results  produced  by  small-amplitude,  first-order  cnoidal,  and 
second-order  cnoidal  wave  theories.  In  Section  A. 3,  results  from  small- 
amplitude  and  second-order  cnoidal  wave  theories  for  a  spherical  shoal 
and  an  idealized  trench  are  compared.  A  discussion  of  the  computer  time 
required  for  second-order  cnoidal  vs.  small-amplitude  wave  modeling  is 
given  in  Section  6.4. 

For  simulations  using  either  first-order  or  second-order  cnoidal 
wave  theory,  the  numerical  model  developed  in  this  report  was  used.  For 
simulations  using  small-amplitude  wave  theory,  a  model  employing 
numerical  methods  similar  to  the  cnoidal  model  was  developed.  The  wave 
angle  was  determined  from  the  irrotionality  of  wave  number,  v  x  k  =  0  , 
and  wave  height  was  determined  by  using  conservation  of  energy  flux, 
v  •  F  =  0.  Equations  6.1  and  6.2  were  used  to  calculate  k  and  F  for 
small -amplitude  wave  theory.  (The  subscript  SA  denotes  small  amplitude.) 


It  is  customary  to  identify  waves  by  the  deepwater  values  of  LQ  and  H0, 
as  determined  by  small -amplitude  wave  theory.  Since  cnoidal  wave  theory 
is  not  valid  in  deep  water,  boundary  conditions  are  needed  on  the  sea¬ 
ward  boundary  of  the  cnoidal  wave  model  if  the  simulations  are  to  be 
characterized  by  deepwater  wave  parameters.  In  Chapter  I  it  was 
mentioned  that  a  future  version  of  the  second-order  cnoidal  wave  model 
would  be  coupled  with  a  model  using  third-order  Stokes  wave  theory  for 
use  in  deeper  water.  The  coupling  of  these  two  models  will  require 
further  research.  For  the  results  presented  in  this  chapter,  the  small- 
amplitude  wave  model  described  above  was  used  to  transform  the  deepwater 
wave  inputs  (HQ  and  ao)  to  the  seaward  boundary  of  the  cnoidal  wave 
model.  The  seaward  boundary  of  the  cnoidal  model  was  placed  at  the  depth 
where  the  value  of  the  Ursell  parameter  (U  =  HL2/d3)  was  equal  to  15,  as 
calculated  from  small-amplitude  theory.  Thus,  wave  heights  and  wave 
angles  were  identical  at  the  "connection  point"  for  all  three  theories 
(small-amplitude,  first-order  cnoidal,  and  second-order  cnoidal)  if  the 
same  deepwater  wave  condition  is  simulated. 

In  all  the  simulations  conducted  in  this  chapter,  a  square  grid 
spacing  of  5  m  was  used.  The  acceleration  of  gravity  was  set  equal  to 
9.806  m/s2,  and  the  density  of  sea  water  was  set  at  1026  kg/m^. 

6.1  Shoaling  Over  a  Plane  Bottom 

In  Section  6.1.1  the  shoaling  of  small -amplitude,  first-order 
cnoidal,  and  second-order  cnoidal  waves  is  compared.  Numerical  results 
from  the  three  wave  theories  are  compared  to  laboratory  wave  flume 
results  in  Section  6.1.2.  Plots  for  multiple  combinations  of  wave  period 
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and  deepwater  wave  steepness  comparing  second-order  cnoidal  wave  theory 
to  small-amplitude  wave  theory  are  presented  in  Section  6.1.3. 

6.1.1  Comparison  of  wave  shoaling  among  small-amplitude,  first- 
order  cnoidal,  and  second-order  cnoidal  waves 

Small-amplitude,  first-order  cnoidal,  and  second-order  cnoidal 
wave  theories  are  compared  for  the  calculation  of  shoaling  in  Figures 

6.1  and  6.2.  Results  are  presented  for  two  deepwater  wave  steepnesses 
(Hq/L0  equal  to  0.005  and  0.020)  at  each  of  two  wave  periods  (8  and 

12  s).  The  ordinate  of  these  plots  is  nondimens ional  wave  height 
(H/Hq),  and  the  abscissa  is  depth  in  meters.  The  curves  stop  where  the 
breaking  criterion  H/D  =  0.8  is  reached. 

Several  features  can  be  seen  in  Figures  6.1  and  6.2.  Second-order 
cnoidal  wave  theory  (cn  II)  always  produces  a  shoaling  rate  intermediate 
between  the  rates  predicted  by  the  two  other  theories.  This  was  found 
to  be  true  in  all  cases  simulated.  The  highest  shoaling  rate  is  always 
given  by  first-order  cnoidal  wave  theory  (cn  I),  and  the  lowest  rate  by 
small -amplitude  wave  theory  (SA).  For  a  given  wave  period,  the  cn  II 
curve  is  closest  to  the  cn  I  curve  at  the  smaller  deepwater  wave 
steepness  (Figures  6.1a  and  6.2a)  and  closest  to  the  SA  curve  at  the 
greater  deepwater  wave  steepness  (Figures  6.1b  and  6.2b). 

As  described  in  Chapter  V,  there  are  three  steps  in  the  calcu¬ 
lation  of  wave  height  using  cnoidal  wave  theory: 

a.  The  energy  flux,  F,  is  calculated  for  the  seaward  boundary 

based  on  the  primary  variables  H,  T,  and  D  (Equations  5.15  and 
5.16).  Equations  5.15  and  5.16  are  repeated  as 
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Figure  6.1  Comparison  of  wave  shoaling:  SA,  on  I,  and 
cn  II  waves,  T  =  8  s. 
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=  H2pg‘)JgD^F0  ♦  eF^ 


Recall  that  F  for  first-order  cnoidal  wave  theory  is 
calculated  from  Equation  6.3  if  eF^  is  set  equal  to  zero. 

b.  For  the  case  of  pure  shoaling,  i.e.,  aQ  =  0  ,  the  value  of 

the  energy  flux  remains  constant  at  each  grid  node  as  the  wave 

propagates  towards  the  point  of  wave  breaking. 

c.  The  energy  flux  has  two  components,  Hc  and  F  .  As  F  which 
is  a  function  of  Ursell  parameter  and  depth,  changes  value, 
the  wave  height  is  calculated  by  imposing  the  conservation 

of  energy  flux  (Equation  5.17).  Equation  5.17  is  repeated  as 


H 


1/2 


(6.4) 


Three  terms  (  Vi°*’  Fq,  and  eFq)  in  the  equation  for  F* 

(Equation  6.3)  change  in  value  as  the  wave  shoals.  The  terms -^gD”  and  Fq 
decrease  during  shoaling,  exerting  a  tendency  to  produce  larger  wave 
heights.  However,  the  third  terra,  eF^,  which  occurs  in  the  calculation 
of  second-order  energy  flux,  increases  as  the  wave  shoals,  negating  some 
of  the  effect  of  the  first  two  terms.  The  value  of  e  is  directly 
affected  by  the  value  of  the  deepwater  wave  steepness.  The  greater  the 
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value  of  Hq/L0,  the  greater  the  value  of  e  throughout  the  domain 
of  the  simulation,  and  the  greater  the  effect  of  the  eF^  term  in  the 

n 

calculation  of  F  .  Thus,  the  shoaling  curve  for  cn  II  is  "pushed 
downward,"  away  from  the  cn  I  shoaling  curve,  with  the  difference 
between  the  two  cnoidal  theories  increasing  at  larger  deepwater  wave 
steepnesses. 

Multiple  wave  period  and  deepwater  wave  steepness  combinations 
were  simulated  for  the  pure  shoaling  case  (aQ  =  0}  over  a  plane  bottom 
slope  of  1:50.  For  each  of  five  wave  periods  (6,  8,  10,  12,  and  14  s), 
six  deepwater  wave  steepnesses  (H0/LQ  =  0.005,  0.010,  0.015,  0.020, 
0.025,  0.030)  were  simulated.  The  ensemble  of  simulations  was  repeated 
employing  small -amplitude,  first-order  cnoidal,  and  second-order  cnoidal 
wave  theories.  Table  6.1  contains  the  wave  height  at  breaking  (Hb)  and 
the  location  of  wave  breaking  relative  to  the  shoreline  (xb)  calculated 
by  each  of  the  three  wave  theories.  Table  6.2  contains  of  the  maximum 
and  minimum  percent  differences  among  the  three  theories  calculated  from 
the  results  shown  in  Table  6.1.  Also  in  Table  6.2  are  the  values  of 
Hq/L0  corresponding  to  the  maximum  and  minimum  percent  differences.  The 
percent  difference  between  any  two  wave  theories  is  a  function  of  deep¬ 
water  wave  steepness  but  is  not  a  function  of  wave  period  for  a  given 
deepwater  wave  steepness. 

In  Table  6.2  the  largest  percent  differences  between  SA  and  cn  II 
occur  at  HQ/L0  =  0.005,  with  Hb  and  xb,  determined  from  cn  II,  22% 
greater  than  the  SA  values.  This  difference  diminishes  to  3%  at  VLo 
=  0.030.  The  largest  percent  difference  between  cn  I  and  SA  is  also  at 
H0/Lq  =  0.005,  with  the  cn  I  values  approximately  37%  higher 
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Table  6.1.  Shoaling  over  a  plane  bottom  (1:50):  Height 
and  position  of  waves  at  breaking  (Hb  =  0.8D) 
as  determined  by  SA,  cn  I,  and  cn  11° 


Hb  xb 


T 

Ho/Lo 

Hd 

(a) 

(a) 

(s) 

(a) 

SA 

cn  I 

cn  II 

SA 

cn  I 

cn  II 

0.005 

0.28 

0.41 

0.56 

0.50 

26 

35 

32 

0.010 

0.56 

0.72 

0.92 

0.82 

45 

SB 

52 

A 

0.015 

0.84 

1.01 

1.25 

1.11 

63 

78 

69 

0.020 

1.12 

1.28 

1.55 

1.37 

80 

97 

86 

0.025 

1.40 

1.54 

1.84 

1.62 

96 

115 

101 

0.030 

1.69 

1.79 

2.11 

1.85 

112 

132 

116 

0.005 

0.50 

0.73 

1.00 

0.89 

46 

63 

56 

0.010 

1.00 

1.29 

1.64 

1.47 

80 

103 

92 

8 

0.015 

1.50 

1.79 

2.22 

1.97 

112 

139 

123 

0.020 

2.00 

2.27 

2.75 

2.43 

142 

172 

152 

0.025 

2.50 

2.74 

3.26 

2.87 

171 

204 

180 

0.030 

3.00 

3.19 

3.75 

3.29 

199 

234 

206 

0.005 

0.78 

1.14 

1.56 

1.39 

72 

98 

87 

0.010 

1.56 

2.01 

2.57 

2.29 

126 

161 

143 

10 

0.015 

2.34 

2.80 

3.46 

3.08 

175 

217 

192 

0.020 

3.12 

3.55 

4.29 

3.80 

222 

269 

238 

0.025 

3.90 

4.28 

5.09 

4.49 

267 

318 

281 

0.030 

4.68 

4.98 

5.86 

5.14 

312 

366 

321 

0.005 

1.12 

1.65 

2.25 

2.00 

103 

140 

125 

0.010 

2.25 

2.89 

3.70 

3.30 

181 

231 

206 

12 

0.015 

3.37 

4.03 

4.99 

4.43 

252 

312 

277 

0.020 

4.49 

5.11 

6.19 

5.48 

320 

387 

342 

0.025 

5.62 

6.16 

7.33 

6.46 

385 

458 

404 

0.030 

6.74 

7.18 

8.43 

7.40 

449 

527 

463 

0.005 

1.53 

2.24 

3.06 

2.73 

140 

191 

170 

0.010 

3.06 

3.94 

5.04 

4.49 

246 

315 

281 

14 

0.015 

4.59 

5.49 

6.79 

6.03 

343 

424 

377 

0.020 

6.12 

6.96 

8.42 

7.45 

435 

526 

466 

0.025 

7.65 

8.38 

9.97 

8.80 

524 

623 

550 

0.030 

9.18 

9.77 

11.48 

10.07 

611 

718 

630 

i 
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Table  6.2.  Shoaling  over  a  plane  bottom  (1:50):  Maximum  and 
minimum  percent  differences  of  wave  height  at 
breaking  among  SA,  cn  I,  and  cn  II 


HAX 

Breaking 

at 

Height 

I1IN 

at 

XDIFF 

fc>/Lo 

IDIFF 

Ho/Lo 

cn  I  :  SA 

37 

0.005 

17 

0.030 

cn  II  :  SA 

22 

0.005 

3 

0.030 

cn  I  :  cn  II 

14 

0.030 

12 

0.010 

than  the  corresponding  SA  values.  This  difference  diminishes  to  18%  at 
Hq/L0  =  0.030.  The  difference  between  Hb  as  determined  from  cn  I  and  cn 
II  shows  the  opposite  trend;  the  largest  percent  difference  (as  was 
noted  in  Figures  6.1  and  6.2)  occurs  at  HQ/L0  =  0.030,  with  cn  I  values 
14%  higher  than  cn  II  values.  This  difference  is  relatively  constant 
over  the  range  of  deepwater  wave  steepness  considered,  reducing  only  to 
12%  at  H0/Lq  =  0.005. 

6.1.2  Comparison  of  Theoretical  and  Experimental  Shoaling  Rates 

Numerical  model  results  for  shoaling  over  a  plane  beach  were 
compared  to  laboratory  experimental  data  from  Buhr  Hansen  and  Svendsen 
(1979).  In  this  set  of  experiments,  shoaling  on  a  plane  bottom  (1:34.26) 
for  waves  of  various  deepwater  wave  steepnesses  was  investigated.  The 
experimenters  suppressed  the  contaminating  free  second  harmonic 
disturbance,  which  usually  results  from  the  difference  between  the 
motion  of  the  wave  generator  and  the  water  particle  velocities  of  waves, 
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by  driving  the  generator  with  a  signal  that  contained  a  component 
designed  to  cancel  the  second  harmonic.  A  profile  view  of  the  wave 

flume  used  in  the  experiment  is  shown  in  Figure  6.3. 

I 

'  The  results  of  17  tests  were  published.  Four  of  these  tests  were 

chosen  for  comparison  with  theoretical  shoaling  predicted  by  the 
numerical  model.  The  four  tests  (Table  6.3)  were  chosen  to  cover  a 
range  of  deepwater  wave  steepnesses  (0.004  to  0.025)  for  tests 
applicable  to  cnoidal  wave  theory  (  U  >10).  Numerical  results  were 
calculated  starting  at  the  beginning  of  the  sloping  bottom  of  the  wave 
flume.  The  wave  height  for  input  to  the  numerical  models  was  the  wave 
height  measured  at  this  point.  Note  that  the  deepwater  wave  steepnesses 
shown  in  Table  6.3  differ  from  those  reported  by  the  experimenters. 

The  deepwater  wave  steepnesses  in  Table  6.3  were  calculated  using 

I 

I 

t 


0.36m 


Figure  6.3.  The  wave  flume  used  in  the  Buhr  Hansen  and 
Svendsen  (1979)  shoaling  experiments 

I 
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small-amplitude  wave  theory  and  the  wave  height  measured  at  the 
beginning  of  the  sloping  portion  of  the  wave  flume. 


Table  6.3.  Experimental  data  of  wave  shoaling 


Test* 

H 

(at  start  of  slope) 

(m) 

H  /L 
o  o 

VDb 

(data) 

H 

0.097 

0.025 

0.937 

K 

0.080 

0.021 

0.861 

N 

0.065 

0.011 

0.982 

Q 

0.039 

0.004 

1.001 

*  Letter  identifying  the  test  corresponds  to  the  notation 
of  Buhr  Hansen  and  Svendsen  (1979). 

The  results  from  the  numerical  simulations  (SA,  cn  I,  and  cn  II) 
are  compared  to  the  experimental  results  in  Figures  6.4  to  6.7.  The 
numerical  simulations  were  continued  until  the  wave  height  to  depth 
ratio  at  breaking  measured  in  the  flume  experiments  (Table  6.3)  was 
reached . 

The  same  relation  among  the  three  theoretical  shoaling  curves  that 
was  noted  in  Figures  6.1  and  6.2  is  seen  in  Figures  6.4  to  6.7  .  The  cn 
II  curve  is  situated  between  the  higher-lying  cn  I  curve  and  the  lower- 
lying  SA  curve.  The  cn  II  curve  is  closer  to  the  cn  I  curve  at  smaller 
deepwater  wave  steepnesses  and  closer  to  the  SA  curve  at  larger 
steepnesses.  For  all  four  cases,  the  cn  II  curve  more  accurately 
describes  the  measurements  than  the  SA  curve.  For  the  two  larger 
deepwater  wave  steepnesses  (Figures  6.4  and  6.5),  the  SA  curve  is  a 
better  match  to  the  data  than  the  cn  I  curve,  except  in  the  region  Just 
prior  to  wave  breaking,  where  the  SA  curve  grossly  underpredicts  the 
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Figure  6.5  Comparison  with  flume  data  -  Test  K 
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Figure  6.6  Comparison  with  flume  data  -  Test  N 
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Figure  6.7  Comparison  with  flume  data  -  Test  Q 
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trend  of  the  data.  The  shoaling  rate  is  overpredicted  by  cn  I  in  deeper 
water  for  the  two  larger  deepwater  wave  steepnesses  (Figures  6.4  and 
6.5),  and  although  the  cn  I  curve  approaches  the  experimental  data  at 
the  breaking  point  in  these  two  figures,  the  shape  of  the  experimental 
curve  is  not  duplicated.  The  cn  II  curve  is  an  acceptable  match  to  the 
experimental  data  in  all  four  cases  except  near  the  breaking  point, 
where  the  cn  II  curve  underpredicts  the  trend  of  the  data,  especially 
for  the  larger  deepwater  wave  steepnesses  (Figures  6.4  and  6.5). 

The  failure  of  second-order  cnoidal  wave  theory  to  accurately 
match  experimental  data  near  to  wave  breaking,  although  somewhat 
disappointing,  should  not  be  completely  unexpected.  As  discussed  in 
Section  6.1.1  the  term  eF^  exerts  a  downward  influence  on  the  shoaling 
curve  which  increases  as  e  =  H/D  increases.  A  perturbation  solution 
which  assumes  that  e  is  small  can  not  be  expected  to  perform  well  if 
the  value  of  e  approaches  1.0.  These  results  suggest  that  third-order 
cnoidal  wave  theory  may  provide  better  results  near  wave  breaking  and 
for  steeper  waves. 


6.1.3  The  shoaling  of  second-order  cnoidal  waves  over  a  plane  bottom 


Shoaling  calculated  by  second-order  cnoidal  wave  theory  is 
compared  to  shoaling  calculated  by  small-amplitude  wave  theory  in 
Figures  6.8  to  6.12.  The  bottom  slope  used  for  all  cases  was  1:50.  Six 
deepwater  wave  steepnesses  (HQ/L0  =  0.005,  0.010,  0.015,  0.020,  0.025, 
and  0.030)  were  simulated  for  each  of  five  wave  periods  (T  =  6,  8,  10, 
12,  and  14  s).  The  curve  for  H0/L0  =  0.030  was  not  plotted  since  it  was 
indistinguishable  from  the  curve  for  HQ/L0  =  0.025  .  The  ordinate  in 
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Figure  6.8  Shoaling  of  second-order  cnoidal  ard  small -amplitude 
waves  over  a  plane  bottom  (1:50)  T  =  6  s. 
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Figure  6.9.  Shoaling  of  second-order  cnoidal  and  small -amplitude 
waves  over  a  plane  bottom  (1:50),  T  =  8  s. 
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Figure  6.10.  Shoaling  of  second-order  cnoidal  and  small-amplitude 
waves  over  a  plane  bottom  (1:50),  T  =  10  s. 
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Figure  6.11.  Shoaling  of  second-order  cnoidal  and  small -amplitude 
waves  over  a  plane  bottom  (1:50),  T  =  12  s. 


Figure  6.12.  Shoaling  of  second-order  cnoidal  and  small-amplitude 
waves  over  a  plane  bottom  (1:50),  T  =  14  s. 

these  plots  is  nondimens ional  wave  height  (H/HQ).  The  abscissa  is  water 
depth  in  meters,  and  the  scale  of  the  abscissa  varies  from  plot  to 
plot.  The  curves  for  the  five  deepwater  wave  steepnesses  on  each  plot 
are  halted  at  reaching  the  breaking  criterion  of  H/D  =0.8  .  The  small- 
amplitude  curve  shown  in  Figures  6.8  to  6.12  is  not  halted  at  a  breaking 
point.  Small -amplitude  shoaling  is  not  dependent  on  wave  steepness,  for 
a  given  period;  therefore,  only  one  SA  curve  need  be  presented.  The  SA 
curve  is  continued  past  all  the  cn  II  curves  for  ease  in  comparison. 

6.2  Refraction  Over  a  Plane  Bottom 

Refraction  over  a  plane  bottom  of  small-amplitude,  first-order 
cnoidal,  and  second-order  cnoidal  waves  is  investigated  in  Section 
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6.2.1.  Plots  comparing  small-amplitude  and  second-order  cnoidal  wave 
refraction  are  presented  in  Section  6.2.2  . 

6.2.1  Comparison  of  wave  refraction  among  small-amplitude,  first- 
order  cnoidal,  and  second-order  cnoidal  waves 

The  refraction  of  small-amplitude,  first-order  cnoidal,  and 
second-order  cnoidal  waves  is  compared  in  Figures  6.13  and  6.14  for  wave 
periods  of  6  and  14  s,  respectively.  For  each  figure,  the  deepwater 
wave  angle  aQ  =  60°  and  the  bottom  slope  is  1:50.  In  Figures  6.13a  and 
6.14a,  Hq/L0  =  0.005  and  in  Figures  6.13b  and  6.14b,  HQ/L0  =  0.020. 

Note  that  the  curves  plotted  in  these  figures  are  not  traces  of  the  wave 
rays.  The  ordinate  is  nondimens ional  wave  angle  a/aQ,  and  the  abscissa 
is  water  depth  in  meters.  The  two  cnoidal  simulations  began  with  the 
wave  height,  wave  angle,  and  depth  where  U  =  15,  as  calculated  using 
small-amplitude  wave  theory.  The  plots  stop  where  the  breaking 
criterion  H/D  =  0.8  is  reached. 

As  in  the  shoaling  results,  the  cn  II  curve  lies  between  the  cn  I 
curve  and  the  SA  curve.  Less  refraction  is  predicted  by  cn  I  and 
greater  refraction  is  predicted  by  SA. 

The  cn  I  curves  show  an  interesting  characteristic  Just  seaward 
of  the  wave  breaking  point.  The  change  in  wave  angle  "stalls,"  and  the 
angle  increases  very  slightly  with  decreasing  depth  for  the  HQ/L0 
=  0.005  cases.  Referring  to  Snell's  law  (Equation  4.8),  this  behavior 
requires  that,  near  to  wave  breaking,  the  wavelength  remain  constant 
with  decreasing  depth  ,  and  in  the  case  of  increasing  angle  with 
decreasing  depth,  the  wavelength  must  increase  with  decreasing  depth. 
Svendsen  (1974)  and  Wang  and  Le  Mehaute  (1980)  provide  plots  of 
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Figure  6.13.  Comparison  of  wave  refraction:  SA,  cn  I,  and 
cn  II  waves,  T  =  6  s,  aQ  =  60° . 
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Figure  6.14.  Comparison  of  wave  refraction:  SA,  cn  I,  and 
cn  II  waves,  T  =  14  s,  aQ  =  60°. 


L/LQ  vs  D/Lq  which  show  that  the  wavelength  calculated  using  first- 
order  cnoidal  wave  theory  begins  increasing  shoreward  of  approximately 
H/D  =  0.78  for  a  plane  bottom.  This  tendency  may  indicate  a  limit  of 
applicability  for  first-order  cnoidal  wave  theory,  since  the  increase  of 
wavelength  with  decreasing  depth  is  counter  to  intuition,  and  no  known 
evidence  exists  to  support  this  phenomenon. 

As  is  discussed  in  Chapter  IV,  the  calculation  of  the  angle  of 
wave  propagation  is  important  for  the  determination  of  the  wave-induced 
longshore  current  and  resulting  wave-induced  longshore  sediment  trans¬ 
port.  Therefore,  the  unverified  behavior  of  refraction  calculations 
near  the  wave  breaking  point  obtained  by  using  first-order  cnoidal  wave 
theory  is  a  strong  incentive  to  use  second-order  rather  than  first-order 
cnoidal  theory  for  refraction  calculations. 

The  refraction  of  shallow  water  waves  was  investigated  in  a  manner 
similar  to  the  investigation  of  shoaling  in  Section  6.1.3  .  A  bottom 
slope  of  1:50  was  used  in  all  the  simulations.  Five  wave  periods  (6,  8, 
10,  12,  and  14  s)  combined  with  each  of  six  deepwater  wave  steepnesses 
(0.005,  0.010,  0.015,  0.020,  0.025,  and  0.030)  were  simulated.  This 
ensemble  of  thirty  simulations  was  repeated  for  each  of  three  deepwater 
wave  angles  (30°,  45°,  and  60°).  (See  Figure  5.1  for  a  definition  of 
wave  angle.)  The  resulting  ninety  simulations  were  repeated  using 
small -amptitude,  first-order  cnoidal,  and  second-order  cnoidal  wave 
theories,  for  a  grand  total  of  270  numerical  simulations.  The  seaward 
boundary  conditions  for  the  simulations  employing  either  first  or 
second-order  cnoidal  wave  theory  were  obtained  from  the  simulations 
using  small-amplitude  wave  theory  at  U  =  15.  The  results  for  deepwater 


wave  angles  30°,  45°,  and  60°  are  presented  in  Tables  6.4,  6.5,  and  6.6, 


respectively.  These  tables  list  the  values  of  the  wave  height  at 


breaking,  the  wave  angle  at  breaking,  and  the  location  of  the  wave 


breaking  point  relative  to  the  shoreline. 


Table  6.7  lists  the  maximum  and  minimum  percent  differences  among 


the  three  wave  theories  in  the  calculated  quantities  of  Hb  and  <*b.  The 


maximum  percent  differences  for  occur  at  the  smallest  deepwater  wave 


steepness  and  the  minimum  percent  differences  at  the  largest  deepwater 


wave  steepness  in  a  comparison  of  SA  against  either  cn  I  or  cn  II.  The 


percent  difference  between  cn  I  and  cn  II  wave  angles  at  breaking  is 


almost  constant  with  deepwater  wave  steepness.  Recalling  the  results 


for  shoaling  presented  in  Table  6.2,  Table  6.7  shows  the  same  trends  for 


differences  in  the  calculation  of  Hb  by  the  three  wave  theories. 


It  is  interesting  to  compare  the  magnitudes  in  the  differences  in 


Hfa  recorded  in  Tables  6.2  and  6.7  .  As  discussed  in  Chapter  V,  the 


numerical  calculation  of  refraction  implicitly  includes  the  effects  of 


shoaling.  The  percent  differences  at  aQ  =  30°  are  very  similar  to  the 


results  for  shoaling  (a  =0°).  The  results  at  a  =  45°  and  a  =  60c 

o  o  o 


show  only  a  small  increase  in  the  percent  differences  among  the  three 


theories  in  the  calculation  of  Hb.  It  must  be  concluded  that  the  per¬ 


cent  difference  in  wave  height  at  breaking  among  the  three  wave  theories 


is  due  primarily  to  differences  in  the  calculation  of  shoaling. 


6.2.2  The  refraction  of  second-order  cnoidal  waves  over  a  plane  bottom 


Figures  6.15  to  6.29  are  plots  (nondimensional  wave  angle  vs 


depth)  for  the  second-order  cnoidal  refraction  simulations.  In  a 
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Table  6.4.  Refraction  over  a  plane  bottom  (1:50):  Height,  angle, 
and  position  of  waves  at  breaking,  a  =  30° 


a/lo  Ho 
(■) 


at 

<deg> 

SA  cn  I  cn  II 

6.7  9.9  8.7 

8.8  12.4  10.8 

10.3  14.0  12.3 

11.6  15.3  13.5 

12.6  16.4  14.4 

13.6  17.3  15.2 


xb 

(■) 

SA  cn  I  cn  I! 

25  34  30 

43  55  49 

60  75  66 

76  93  82 

92  110  97 

107  127  111 


43  60  s: 

76  99  88 

106  133  118 

135  165  146 

163  196  172 

190  226  197 

68  93  83 

119  154  137 

166  208  184 

211  258  228 

254  306  269 

297  353  308 

97  134  120 

171  222  197 

239  299  265 

304  372  326 

366  441  387 

427  506  444 


183  16: 
302  26£ 
407  361 
506  444 
600  52: 
691  60i 


Table  6.5.  Refraction  over  a  plane  bottom  (1:50):  Height,  angle,  and 
position  of  waves  at  breaking,  aQ  =  45° 


at 

xb 

Ho/Lo 

Ho 

(!) 

(deg) 

(■) 

(■) 

SA 

cn  I 

cn  11 

SA 

cn  I 

cn  11 

SA 

cn  I 

cn  I! 

0.005 

0.28 

0.36 

0.50 

0.45 

9.1 

13.7 

12.0 

23 

32 

28 

0.010 

0.56 

0.63 

0.83 

0.74 

12.0 

17.1 

14.9 

40 

52 

46 

0.015 

0.84 

0.89 

1.12 

0.99 

14.1 

19.5 

17.0 

55 

70 

62 

0.020 

1.12 

1.13 

1.40 

1.23 

15.9 

21.4 

18.7 

70 

88 

77 

0.025 

1.40 

1.36 

1.66 

1.45 

17.4 

22.9 

20.0 

as 

204 

91 

0.030 

1.69 

1.59 

1.91 

1.67 

18.7 

24.2 

21.2 

99 

120 

104 

0.005 

0.50 

0.64 

0.89 

0.79 

9.0 

13.6 

11.9 

40 

56 

50 

0.010 

1.00 

1.13 

1.48 

1.31 

12.0 

17.1 

15.0 

70 

92 

82 

0.015 

1.50 

1.57 

1.99 

1.76 

14.1 

19.4 

17.0 

98 

125 

110 

0.020 

2.00 

2.00 

2.48 

2.18 

15.9 

21.3 

18.6 

125 

155 

137 

0.025 

2.50 

2.42 

2.95 

2.58 

17.4 

22.9 

20.0 

151 

184 

161 

0.030 

3.00 

2.82 

3.40 

2.96 

18.7 

24.3 

21.2 

176 

213 

185 

0.005 

0.78 

1.00 

1.39 

1.24 

9.0 

13.6 

11.9 

63 

87 

78 

0.010 

1.56 

1.76 

2.31 

2.04 

12.0 

17.1 

14.9 

110 

144 

128 

0.015 

2.34 

2.46 

3.12 

2.75 

14.1 

19.5 

17.0 

154 

195 

172 

0.020 

3.12 

3.13 

3.88 

3.41 

15.8 

21.4 

18.7 

195 

242 

213 

0.025 

3.90 

3.77 

4.61 

4.03 

17.4 

22.9 

20.0 

236 

288 

252 

0.030 

4.68 

4.41 

5.31 

4.63 

18.7 

24.3 

21.2 

276 

332 

289 

0.005 

1.12 

1.44 

2.01 

1.79 

9.0 

13.7 

12.0 

90 

126 

112 

0.010 

2.25 

2.53 

3.32 

2.94 

12.0 

17.1 

14.9 

158 

208 

184 

0.015 

3.37 

3.54 

4.49 

3.97 

14.1 

19.5 

17.0 

221 

281 

248 

0.020 

4.49 

4.50 

3.59 

4.91 

15.8 

21.3 

18.6 

281 

349 

307 

0.025 

5.62 

5.43 

6.63 

5.81 

17.4 

22.9 

20.0 

340 

415 

363 

0.030 

6.74 

6.35 

7.65 

6.67 

18.7 

24.3 

21.2 

397 

478 

417 

0.005 

1.53 

1.96 

2.73 

2.43 

9.0 

13.7 

12.0 

122 

171 

152 

0.010 

3.06 

3.45 

4.52 

4.01 

12.0 

17.1 

14.9 

216 

282 

250 

0.015 

4.39 

4.82 

6.11 

5.40 

14.1 

19.5 

17.0 

301 

382 

337 

0.020 

6.12 

6.13 

7.60 

6.68 

15.8 

21.3 

18.6 

383 

475 

418 

0.025 

7.65 

7.40 

9.03 

7.90 

17.4 

22.9 

20.0 

462 

564 

494 

0.030 

9.18 

8.64 

10.41 

9.07 

18.7 

24.3 

21.2 

540 

651 

567 

97 


Table  6.6.  Refraction  over  a  plane  bottom  (1:50):  Height,  angle,  and 
position  of  waves  at  breaking,  aQ  =  60° 


It) 

ab 

xb 

T 

Ho/lo 

Ho 

(■) 

(deg) 

(■) 

(s) 

(■) 

SA 

cn  I 

cn  11 

SA 

cn  I 

cn  II 

SA 

cn  I 

cn  II 

0.003 

0.28 

0.31 

0.44 

0.40 

10.4 

16.0 

13.9 

20 

28 

25 

0.010 

0.56 

0.55 

0.74 

0.65 

13.7 

20.0 

17.5 

35 

46 

41 

6 

0.015 

0.84 

0.77 

1.00 

0.88 

16.2 

22.9 

19.9 

48 

63 

55 

0.020 

1.12 

0.98 

1.24 

1.09 

18.3 

25.2 

21.9 

61 

78 

68 

0.025 

1.40 

1.18 

1.47 

1.29 

20.0 

27.0 

23.5 

74 

92 

81 

0.030 

1.69 

1.38 

1.70 

1.48 

21.6 

28.7 

25.0 

87 

106 

93 

0.005 

0.50 

0.56 

0.79 

0.70 

10.4 

15.9 

13.9 

35 

50 

44 

0.010 

1.00 

0.98 

1.31 

1.16 

13.7 

20.0 

17.5 

61 

82 

73 

8 

0.015 

1.50 

1.37 

1.77 

1.56 

16.2 

22.8 

19.9 

86 

til 

98 

0.020 

2.00 

1.74 

2.21 

1.94 

18.3 

25.2 

21.9 

109 

138 

121 

0.025 

2.50 

2.11 

2.62 

2.29 

20.0 

27.1 

23.5 

132 

164 

143 

0.030 

3.00 

2.46 

3.03 

2.63 

21.6 

28.8 

25.0 

154 

189 

164 

0.003 

0.78 

0.87 

1.24 

1.10 

10.4 

15.9 

13.9 

54 

78 

69 

0.010 

1.56 

1.53 

2.05 

1.81 

13.7 

20.0 

17.4 

96 

128 

113 

10 

0.015 

2.34 

2.14 

2.77 

2.44 

16.2 

22.8 

19.9 

134 

173 

153 

0.020 

3.12 

2.72 

3.45 

3.03 

18.3 

25.2 

21.9 

170 

215 

189 

0.025 

3.90 

3.29 

4.10 

3.58 

20.0 

27.1 

23.6 

206 

256 

224 

0.030 

4.68 

3.84 

4.73 

4.11 

21.6 

28.8 

25.0 

240 

296 

257 

0.005 

1.12 

1.25 

1.78 

1.58 

10.4 

15.9 

13.9 

78 

112 

99 

0.010 

2.25 

2.21 

2.95 

2.61 

13.7 

20.0 

17.5 

138 

184 

163 

12 

0.015 

3.37 

3.06 

3.99 

3.52 

16.2 

22.8 

19.9 

193 

249 

220 

0.020 

4.49 

3.92 

4.96 

4.36 

18.3 

25.2 

21.9 

245 

310 

272 

0.025 

5.62 

4.74 

5.90 

5.15 

20.0 

27.1 

23.5 

296 

369 

322 

0.030 

6.74 

5.54 

6.81 

5.92 

21.6 

28.8 

25.0 

346 

426 

370 

0.005 

1.53 

1.70 

2.43 

2.15 

10.4 

15.9 

13.9 

107 

152 

135 

0.010 

3.06 

3.00 

4.01 

3.55 

13.7 

20.0 

17.4 

188 

251 

222 

14 

0.015 

4.59 

4.20 

5.43 

4.79 

16.2 

22.8 

19.9 

262 

339 

299 

0.020 

6.12 

5.34 

6.75 

5.93 

18.3 

25.2 

21.9 

334 

422 

371 

0.025 

7.65 

6.45 

8.03 

7.01 

20.0 

27.1 

23.6 

403 

502 

438 

0.030 

9.18 

7.54 

9.27 

8.06 

21.6 

28.8 

25.0 

471 

579 

504 

98 
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Table  6.7.  Refraction  over  a  plane  bottom  (1:50):  Maximum  and  min¬ 
imum  percent  differences  of  wave  height  and  angle  at  breaking  (Hb 


:  0.8D) 

among  SA, 

cn 

I,  and 

cn  II 

(degl 

Rodels 

Capared 

HM 

IDIFF 

Breaking  Hei$>t 
at  HIN 

Ho/la  IDIFF 

at 

Ho/lo 

RAX 

IDIFF 

Breaking  Angle 
at  RIN 

Ho/lo  IDIFF 

at 

Ho/Lo 

cn  I  :  SA 

38 

0.003 

19 

0.030 

30 

0.003 

26 

0.030 

30 

cn  II  :  SA 

23 

0.003 

4 

0.030 

32 

0.005 

12 

0.030 

cn  I  :  cn  II 

14 

0.030 

10 

0.003 

13 

0.010 

13 

0.030 

cn  I  :  SA 

40 

0.005 

20 

0.030 

52 

0.005 

29 

0.030 

43 

cn  II  :  SA 

23 

0.003 

3 

0.030 

33 

0.005 

13 

0.030 

cn  I  :  cn  II 

13 

0.030 

11 

0.003 

13 

0.010 

14 

0.010 

cn  I  :  SA 

43 

0.005 

23 

0.030 

54 

0.005 

33 

0.030 

60 

cn  II  :  SA 

29 

0.005 

7 

0.030 

34 

0.005 

16 

0.030 

cn  I  :  cn  II 

13 

0.030 

10 

0.003 

13 

0.023 

14 

0.010 

similar  manner  to  the  plots  for  wave  shoaling,  curves  for  five  deepwater 
wave  steepnesses  (0.005,  0.010,  0.015,  0.020,  0.025)  are  plotted  with  a 
small -amplitude  curve  for  each  of  five  wave  periods  (6,  8,  10,  12,  and 
14  s).  Figures  6.15  to  6.19  are  the  plots  for  T  =  6  s  to  T  =  14  s  for  a 
deepwater  wave  angle  of  30°.  Likewise,  Figures  6.20  to  6.24  and  Fig¬ 
ures  6.25  to  6.29  are  plots  for  aQ  =  45°  and  60°,  respectively.  The 
curves  for  the  various  deepwater  wave  steepnesses  stop  at  the  breaking 
criterion  H/D  =  0.8.  As  with  shoaling,  small-amplitude  wave  refraction 
does  not  depend  on  wave  steepness,  so  only  one  small -amplitude  curve  is 
presented  for  each  wave  period,  and  that  curve  was  allowed  to  continue 
beyond  the  breaking  criterion.  Companion  plots  to  these  wave  angle 
plots  (nondimensional  wave  height  vs  depth)  are  included  in  Appendix  D. 

From  this  series  of  figures  it  is  seen  that  second-order  cnoidal 


waves  refract  less  than  do  small-amplitude  waves.  The  rate  of  refrac¬ 
tion  decreases  with  increasing  deepwater  wave  steepness  for  a  given  wave 
period  for  second-order  cnoidal  waves. 

6.3  Shoaling  and  refraction  over  non-plane  bathymetry 

The  results  presented  in  the  beginning  sections  of  this  chapter 
demonstrate  the  properties  of  shoaling  and  refraction  of  shallow  water 
waves  as  predicted  by  small-amplitude,  first-order  cnoidal,  and  second- 
order  cnoidal  wave  theories.  Even  though  the  results  corresponded  to  a 
plane  bathymetry,  the  physical  mechanisms  which  cause  shoaling  and 
refraction  do  not  change  if  the  bathymetry  becomes  non-planar.  The 
differences  in  the  prediction  of  shoaling  and  refraction  among  the  three 
theories  remain  the  same  as  for  the  non-planar  case.  This  is  true  as 
long  as  the  basic  assumptions  in  the  derivation  of  the  wave  theories, 
especically  the  assumption  that  the  waves  are  of  constant  form,  remain 
reasonably  valid.  Each  location  can  be  thought  of  as  locally  having  a 
plane  bottom,  and  the  shoaling  and  refraction  predicted  by  the  three 
theories  being  determined  by  the  local  bottom  slope.  This  argument,  of 
course,  holds  only  if  the  effects  of  shoaling  and  refraction  alone  are 
considered.  As  the  bottom  becomes  increasingly  complex,  other  transfor¬ 
mation  mechanisms  neglected  in  the  present  model,  notably  diffraction, 
become  increasingly  important. 

Therefore,  given  the  arguments  in  the  preceeding  paragraph,  the 
prediction  of  shoaling  and  refraction  for  a  non-planar  bottom  is  a 
problem  in  the  numerical  calculation,  rather  than  a  problem  of  different 
physics.  The  techniques  used  for  wave  angle  and  wave  height 


Figure  6.15.  Refraction  over  a  plane  bottom  (1:50),  aQ  =  30°,  T  =  6  s. 


Figure  6.17.  Refraction  over  a  plane  bottom  (1:50),  a  =  30°, 

T  =  10  s.  ° 


Figure  6.18.  Refraction  over  a  plane  bottom  (1:50),  a  =  30°, 

T  =  12  s.  0 
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Figure  6.19.  Refraction  over  a  plane  bottom  (1:50),  a  =  30°, 

T  =  14  s.  ° 


DEPTH  Cm) 

Figure  6.20.  Refraction  over  a  plane  bottom  (1:50),  a  =  45°, 

T  =  6  s.  ° 
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Figure  6.21.  Refraction  over  a  plane  bottom  (1:50),  a  =  45° 

T  =  8  s.  ° 
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Figure  6.22.  Refraction  over  a  plane  bottom  (1:50) 
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Figure  6.23.  Refraction  over  a  plane  bottom  (1:50),  o  =  45°, 

T  =  12  s.  0 


Figure  6.24.  Refraction  over  a  plane  bottom  (1:50),  a  =  45°, 

T  =  14  s.  ° 


105 


o/o„  o/a 


Figure  6.25.  Refraction  over  a  plane  bottom  (1:50),  a  =  60' 

T  =  6  s.  ° 


DEPTH  (m) 


Figure  6.26.  Refraction  over  a  plane  bottom  (1:50),  =  60 


a/a„  a/a 


Figure  6.27.  Refraction  over  a  plane  bottom  (1:50),  a  =  60°, 

T  =  10  s.  ° 


DEPTH  (m) 

Figure  6.29.  Refraction  over  a  plane  bottom  (1:50),  a  =  60°, 

T  =  14  s.  ° 

calculations,  derived  in  Chapter  IV  and  implemented  for  the  numerical 
model  in  Chapter  V  should  therefore  still  apply.  These  techniques  were 
tested  over  an  idealized  non-planar  bathymetry.  As  a  measure  of  the 
performance  of  the  numerical  model,  the  results  from  the  shoaling  and 
refraction  over  an  idealized  bathymetry  can  be  anticipated,  and  the 
model  results  can  be  compared  to  these  expectations. 

Model  results  for  the  shoaling  and  refraction  of  waves  over  a 
spherical  shoal  and  an  idealized  trench  are  presented  in  Sections  6.3.1 
and  6.3.2,  respectively.  In  both  cases  the  lateral  boundaries  were 
placed  so  that  the  lateral  boundary  conditions  would  not  affect  the 
results  of  the  simulations. 


6.3.1  Shoaling  and  refraction  over  a  spherical  shoal 


A  lens  with  the  shape  of  a  top  section  of  a  sphere  was  placed  on  a 
bottom  with  otherwise  constant  depth  (9*9  m).  The  numerical  grid  was 
centered  about  the  shoal.  The  maximum  decrease  in  water  depth  caused  by 
the  spherical  shoal  was  2  m,  and  the  radius  of  the  shoal  was  150  m. 

Waves  of  two  deepwater  wave  steepnesses,  HQ/L0  =  0.005  and  HQ/L0 
=  0.020,  for  a  wave  period  of  12  s  were  simulated  using  the  second-order 
cnoidal  model.  For  comparison,  a  12-s  small  amplitude  wave  was  also 
simulated.  Note  that  small -amplitude  shoaling  and  refraction  do  not 
depend  upon  wave  steepness;  therefore,  a  single  SA  simulation  can  be 
used  in  comparisons  with  the  two  wave  steepnesses  simulated  by  the 
cnoidal  model.  A  constant  wave  height  and  an  angle  normal  to  the 
shoreline  were  the  inputs  to  the  seaward  boundary.  The  wave  heights  at 
the  seaward  boundary  were  determined  from  the  small -amplitude  shoaling 
simulations  described  in  Section  6.1  .  This  implies  a  plane  sloping 
bottom  from  deep  water  to  the  boundary  depth  of  9-9  m,  where  the  bottom 
becomes  horizontal  except  for  the  spherical  disturbance. 

The  results  from  simulations  employing  second-order  cnoidal  wave 
theory  are  presented  in  Figures  6.30  and  6.31,  for  HQ/L0  equal  to  0.005 
and  0.020,  respectively.  Results  are  presented  for  three  cross  sections 
parallel  to  the  shoreline.  The  cross  sections  are  located  fifty  meters 
seaward  of  the  center  of  the  sphere,  at  the  center  of  sphere,  and  fifty 
meters  shoreward  of  the  center  of  the  sphere.  There  are  three  parts  to 
each  figure.  At  each  of  the  three  cross  sections,  Figures  6.30a  and 
6.31a  show  the  longshore  variation  in  depth,  Figures  6.30b  and  6.31b 
show  the  longshore  variation  in  wave  height,  and  Figures  6.30c  and  6.31c 


show  the  longshore  variation  wave  angle.  Note  that  in  Figures  6.30a  and 
6.31a,  the  longshore  variation  of  depth  labeled  "50  m  before  centerline" 
is  identical  to  that  labeled  "50  m  after  centerline".  The  figures  are 
orientated  as  if  the  observer  were  at  sea  and  facing  the  shore. 

Positive  and  negative  angles  mean  the  wave  is  going  to  the  left  and 
right,  respectively,  as  in  Figure  5.1. 

The  following  results  were  anticipated: 

a.  Calculated  quantities  should  be  symmetrical  about  a  line 
running  onshore-offshore  at  the  center  of  the  grid. 

b.  The  wave  angle  should  remain  equal  to  zero  on  the  center- 
line  and  at  the  lateral  boundaries  of  each  cross  section. 

c.  At  either  lateral  side  of  the  shoal,  the  wave  angle  should 
be  directed  toward  the  middle  of  the  grid. 

d.  The  maximum  wave  angle  should  be  located  near  the  position 
of  maximum  change  in  the  bathymetry  in  the  longshore 
direction. 

e.  The  maximum  wave  height  should  be  located  at  the  middle  of 
the  grid. 

f.  The  minimum  wave  height  should  be  located  close  to  the 
edge  of  the  spherical  disturbance. 

g.  The  trends  in  items  c  through  e  should  increase  as  the 
wave  propagated  shoreward. 

These  anticipated  results  are  confirmed  by  the  results  of  the 
numerical  simulations  shown  in  Figures  6.30  and  6.31.  The  model 
simulates  waves  propagating  over  the  spherical  disturbance  with  no 
detectable  numerical  difficulties.  Due  to  the  absence  of  combined 
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Figure  6.30  Refraction  over  a  spherical  shoal,  HQ/L0  =  0.005,  T 
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Figure  6.31  Refraction  over  a  spherical  shoal,  HQ/L0  =  0,020,  T  = 
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refraction-diffraction  in  the  numerical  model,  results  in  the  zone  of 
convergence  shoreward  of  the  shoal  will  become  increasingly  inaccurate. 
The  model  did  not  "blow-up"  during  these  runs,  but  did  "blow-up"  when  a 
more  severe  shoal  (5  m)  was  tested.  For  artificial  bathymetries  such  as 
the  spherical  shoal,  the  area  behind  the  shoal  where  energy  will  be 
focused  is  an  obvious  area  where  the  limits  of  applicability  of  the 
model  will  be  tested.  If  real  bathymetries  are  modeled,  these  areas  of 
inapplicability  may  not  be  as  apparent.  For  every  simulation,  model 
results  should  be  scrutinized  to  determine  if  regions  exist  where  either 
unrealistically  high  or  unrealistically  low  wave  heights  are  calculated 
due  to  the  absence  of  diffraction.  These  unrealistic  values  may  not 
always  cause  the  simulation  to  "blow-up,"  triggering  a  direct  sign  of 
model  inapplicability. 

Several  interesting  features  can  be  seen  in  the  Figures  6.30 
and  6.31.  There  is  less  angle  change  with  the  steeper  wave  (Fig¬ 
ures  6. 31c).  This  reproduces  the  results  of  refraction  over  a  plane 
bottom  presented  in  Figures  6.15  to  6.29,  where  at  a  given  depth, 
steeper  waves  refracted  less  than  less  steep  waves.  Looking  at  the 
maximum  wave  heights  on  each  of  the  three  curves  of  Figures  6.30b  and 
6.31b,  the  effects  of  shoaling  vs.  the  effects  of  concentration  of 
energy  due  to  refraction  can  be  seen.  Recall  from  the  results  for 
shoaling  over  a  plane  bottom  (Figures  6.8  to  6.12),  that  at  a  given 
depth  the  steeper  wave  experiences  greater  shoaling.  On  the  seaward 
side  of  the  shoal  both  refraction  and  shoaling  cause  an  increase  in  wave 
height.  On  the  shoreward  side  of  the  shoal,  refraction  continues  to 
bring  increasing  amounts  of  energy  towards  the  center,  producing  a 


greater  wave  height.  In  contrast,  the  increasing  depth  as  the  wave 
propagates  down  the  backside  of  the  lens  will  cause  "de-shoal ing" ,  a 
term  coined  here  to  describe  wave  height  change  resulting  from  wave 
propagation  into  deeper  water.  The  concentration  of  energy  towards  the 
center  of  the  shoal  overwhelms  the  reduction  of  wave  height  due  to  this 
"de-shoal ing,"  resulting  in  an  increased  wave  height  on  the  shoreward 
side  of  the  shoal  for  both  wave  steepnesses.  However,  for  the  steeper 
wave  (Figure  6.31b)  with  its  greater  "de-shoaling"  and  smaller 
refraction,  the  plot  labeled  "50  m  after  centerline"  is  much  closer  to 
the  plot  labeled  "along  centerline"  than  the  corresponding  plots  on  the 
less  steep  wave  (Figure  6.30b). 

A  comparison  between  the  results  of  the  small-amplitude  and 
second-order  cnoidal  models  for  the  spherical  shoal  bathymetry  is 
presented  in  Figure  6.32.  This  figure  is  similar  to  the  previous  two 
figures.  Three  plots  -  one  each  for  depth,  wave  height,  and  wave  angle 
are  included.  The  slightly  uneven  nature  of  the  curves  for  wave  height 
results  from  round-off  error  in  the  number  of  digits  retained  for 
plotting.  The  results  are  shown  for  the  cross  section  at  the  center  of 
the  shoal.  Figure  6.32  echoes  the  results  presented  in  Section  6.2;  the 
curve  for  small-amplitude  wave  shows  less  shoaling  but  more  refraction 
than  do  the  curves  for  the  second-order  cnoidal  waves.  As  would  be 
expected,  the  results  for  the  less  steep  second-order  cnoidal  wave  are 
closer  to  the  results  for  the  small -amplitude  wave  than  are  the  results 
from  the  steeper  second-order  cnoidal  wave. 


6.3.2  Shoaling  and  refraction  over  a  trench 

The  trench  bathymetry  was  calculated  using  Equation  5.1  (a  =  0, 
b  =  0.02,  d  =  1.0  .  Unlike  the  spherical  shoal  example,  the  trench  was 
superimposed  on  a  bottom  with  slope  of  1:50.  A  single  trench  was  placed 
along  the  center  of  the  grid  with  the  trench  axis  in  the  onshore- 
offshore  direction.  The  maximum  increase  in  depth  imposed  on  the 
sloping  plane  bottom  by  the  trench  was  2  m  (for  a  total  depth  of  11  m) 
at  the  centerline  of  the  grid  450  m  offshore.  The  effect  of  the  trench 
on  the  sloping  bottom  was  linearly  increased  from  0  m  at  the  seaward 
boundary  (600  m  offshore)  to  the  maximum  and  then  was  linearly  decreased 
to  0  m  at  the  shoreline.  This  configuration  represents  an  idealized  rip 
channel,  a  feature  commonly  found  along  sandy  coastlines.  The  seaward 
boundary  was  placed  at  a  water  depth  of  12  m.  The  position  of  the 
seaward  boundary  allowed  the  comparison,  at  a  range  of  depths,  for 
simulations  with  deepwater  wave  steepnesses  of  0.005  and  0.020  for  a 
wave  period  of  12  s.  The  smaller  wave  was  inside  the  cnoidal  range 
(U  >  10),  and  the  larger  wave  was  outside  the  point  of  wave  breaking. 

Simulations  were  conducted  with  both  the  small-amplitude  and  the 
second-order  cnoidal  models.  A  constant  wave  height  and  a  wave  direc¬ 
tion  normal  to  the  shoreline  (a  =  0°)  were  the  inputs  on  the  seaward 
boundary.  The  value  of  the  wave  height  on  the  boundary  was  determined 
from  the  the  small -amplitude  shoaling  simulations  presented  in  Sec¬ 
tion  6.1  . 

The  results  of  simulations  employing  the  second-order  cnoidal 
model  are  presented  in  Figures  6.33  and  6.34,  for  HQ/L0  =  0.005  and 
H0/L0  =  0.020,  respectively.  As  with  the  Figures  6.30  and  6.31  for  the 


spherical  shoal,  there  are  three  parts  to  each  figure.  Figures  6.33a 
and  6.34a  show  the  longshore  variation  in  depth  at  three  cross  sections 
(400,  350,  and  300  m  offshore  for  HQ/L0  =  0.005  and  500,  450,  and  400  m 
offshore  for  HQ/L0  =  0.020)  parallel  to  the  shoreline.  Note  that  the 
three  cross  sections  shown  for  the  steeper  wave  span  the  region  of  the 
maximum  increase  in  depth  due  to  the  trench,  whereas  the  cross  sections 
shown  for  the  less  steep  wave  are  all  shoreward  of  the  maximum  trench 
location.  Figures  6.33b  and  6.34b  are  plots  of  the  longshore  variation 
in  wave  height  corresponding  to  the  three  cross  sections  for  the  respec 
tive  deepwater  wave  steepness.  Similarly,  Figures  6.33c  and  6.34c  are 
plots  of  the  longshore  variation  of  wave  angle  at  each  of  the  cross 
sections.  As  with  Figures  6.30  and  6.31,  Figures  6.33  and  6.34  are 
orientated  as  if  the  observer  were  at  sea  facing  toward  the  shore. 

The  following  results  were  anticipated: 

a.  Calculated  quantities  should  be  symmetrical  about  the 
centerline  of  the  trench. 

b.  There  should  be  no  change  in  wave  angle  along  the 
centerline  of  the  trench. 

c.  At  either  side  of  the  trench,  the  wave  angle  should  be 
directed  away  from  the  centerline  of  the  trench. 

d.  The  maximum  wave  angle  should  be  located  near  the 
position  of  maximum  change  in  the  bathymetry  in  the 
longshore  direction. 

e.  The  maximum  wave  height  should  be  located  near  the 
lateral  edges  of  the  trench. 

f.  The  minimum  wave  height  would  be  located  at  the 
centerline  of  the  trench. 


g.  The  trends  in  items  c  through  e  would  increase  as  the 
wave  propagated  shoreward. 

These  anticipated  results  are  confirmed  by  the  results  of  the 
numerical  simulations  shown  in  Figures  6.33  and  6.34.  The  model  simu¬ 
lates  waves  propagating  over  the  idealized  trench  with  no  detectable 
numerical  difficulties. 

A  comparison  between  the  results  of  the  small-amplitude  and 
second-order  cnoidal  models  for  the  trench  bathymetry  is  presented  in 
Figure  6.35.  The  results  shown  are  for  the  cross  section  400  m  offshore. 
Figure  6.35  is  similar  to  the  previous  five  figures.  There  are  three 
parts  to  the  figure,  one  plot  each  for  depth,  wave  height,  and  wave 
angle.  Figure  6.35  repeats  the  results  presented  in  Section  6.2  and 
Figure  6.32;  the  curve  for  small-amplitude  wave  shows  less  shoaling  but 
more  refraction  than  does  the  curve  for  the  second-order  cnoidal  wave. 
The  dependence  of  wavelength  on  wave  height  can  be  seen  in  both  the 
height  and  angle  plots  (Figures  6.35b  and  6.35c).  The  refractive  effect 
of  the  trench  on  the  wave  height  and  angle  of  the  steeper  second-order 
cnoidal  wave  extends  further  onto  the  constant  depth  regions  at  the 
lateral  sides  of  the  grid  than  for  either  the  less  steep  cnoidal  wave  or 
the  small -amplitude  wave.  Likewise,  the  refractive  effect  of  the  trench 
on  the  less  steep  cnoidal  wave  spreads  onto  these  regions  to  a  greater 
extent  than  for  the  small-amplitude  wave.  This  tendency  can  be 
explained  by  the  change  in  wavelength  caused  by  changes  in  wave  height. 
The  wavelength  determined  by  small-amplitude  wave  theory  is  constant 
over  the  constant  depth  regions,  since  wavelength  is  a  function  of  only 
depth  and  wave  period  for  small-amplitude  waves.  However,  at  a  constant 
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depth,  the  wavelength  of  a  second-cnoidal  wave  increases  with  wave 
height  (increasing  e).  This  increase  of  wavelength  due  to  increasing 
wave  height  has  the  same  effect  on  refraction  as  the  increase  in  wave¬ 
length  due  to  increasing  depth.  Thus,  the  finite-amplitude  waves  react 
as  if  the  trench  were  laterally  extended.  The  steeper  wave  shows  this 
effect  more  than  the  less  steep  wave  since  the  steeper  wave  has  a 
greater  value  of  e. 

The  change  in  wave  height  and  angle  with  deepwater  wave  steepness 
predicted  by  second-order  wave  theory  but  unaccounted  for  by  small- 
amplitude  wave  theory  (as  shown  in  Figures  6.30  to  6.35)  is  a  convincing 
argument  for  the  use  of  a  finite-amplitude  wave  theory  to  calculate  wave 
transformation. 

6.4  Comparison  of  Computer  time  between  Small-Amplitude  and 
Second-Order  Simulations 

The  presentation  of  results  would  not  be  complete  without  com¬ 
menting  on  the  costs  of  modeling  finite-amplitude  shallow  water  waves 
using  second-order  cnoidal  wave  theory.  None  of  the  models  used  to 
produce  results  for  this  report  were  tuned  for  maximum  efficiency. 

Array  sizes  were  larger  than  necessary,  and  many  extra  variables  were 
calculated  and  stored  for  research  purposes.  Nevertheless,  some 
conclusions  can  be  reached  regarding  the  relative  costs  of  small- 
amplitude  vs.  finite-amplitude  wave  modeling.  The  simulations  were  run 
on  a  VAX  750  mini-computer,  computing  power  which  is  presently  equalled 
or  exceeded  at  most  universities  and  many  consulting  firms.  The  ratio 
of  CPU  time  between  the  second-order  cnoidal  and  small-amplitude 
simulations  was  approximately  3.5  :  1  .  For  a  simulation  on  a  100  by 


100  grid  this  difference  amounted  to  computing  times  of  Just  over  a 
minute  for  the  small -amplitude  model  and  less  than  four  minutes  for  the 
second-order  cnoidal  model.  Interestingly  the  first-order  cnoidal 
simulations  required  slightly  more  computer  time  than  the  second-order 
cnoidal  simulations,  because  the  greater  shoaling  predicted  by  first- 
order  theory  near  to  the  breaking  point  required  more  iterations  to 
converge. 

If  wave  simulation  required  hours  of  computing  time,  than  the  3.5  :  1 
ratio  would  indeed  be  an  obstacle  to  the  use  of  second-order  cnoidal 
theory.  But  an  extra  few  minutes  will  not  significantly  raise  the  cost 
of  either  computer  or  human  time.  The  extra  expense  incurred  seems  well 
worth  the  extra  accuracy  obtainable  with  a  finite-amplitude  wave  theory. 


CHAPTER  VII 


CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  STUDY 

The  goal  of  this  report  was  the  development  of  a  shallow  water 
wave  transformation  model  that  employs  second-order  cnoidal  wave  theory 
and  has  the  capability  to  directly  calculate  at  numerical  grid  points. 
This  model  was  proposed  to  be  an  improvement  over  existing  methods  for 
calculating  the  transformation  of  monocromatic  finite-amplitude  waves  in 
shallow  water.  A  secondary  goal  was  to  advance  the  use  of  cnoidal  wave 
theory  by  deriving  and  demonstrating  methods  for  its  accurate  and  effi¬ 
cient  calculation. 

A  significant  difference  was  shown  to  exist  between  shoaling  and 
refraction  predicted  by  second-order  cnoidal  wave  theory  and  that  pre¬ 
dicted  by  the  more  commonly  used  small-amplitude  wave  theory.  The 
greater  shoaling  and  lesser  refraction  predicted  by  second-order  cnoidal 
wave  theory  have  important  implications  for  the  calculation  of  many 
engineering  quantities,  in  particular,  for  the  calculation  of  wave- 
induced  currents  and  sediment  transport  in  the  nearshore  zone.  Further¬ 
more,  although  not  as  good  as  had  been  expected,  the  results  of  second- 
order  cnoidal  wave  shoaling  were  found  to  better  match  experimental  data 
than  did  the  results  of  either  small-amplitude  or  first-order  cnoidal 
wave  theories.  The  model  performed  excellently  in  modeling  waves  propa¬ 
gating  over  non-plane  bathymetries .  The  method  employed  for  the  calcu¬ 
lation  of  results  directly  at  numerical  grid  points  is  an  obvious 
improvement  over  the  more  cumbersome  ray  tracing  method.  The  techniques 


for  the  calculation  of  second-order  cnoidal  wave  shoaling  and  refraction 
used  in  the  numerical  model  result  in  simulations  which  should  be  well 


within  any  financial  limits  imposed  in  a  study  requiring  wave  modeling. 
Although  the  second-order  cnoidal  model  requires  more  computer  time  than 
a  model  based  on  small -amplitude  theory,  the  differences  in  time  and 
cost  are  inconsequential  compared  to  the  increase  in  accuracy  of  the 
results. 

There  is  one  effort  that  remains  before  this  finite-amplitude 
model  can  be  used  in  a  production  mode;  investigation  into  the  connec¬ 
tion  with  the  companion  deeper  water  model,  which  uses  Stokes  third- 
order  theory,  must  be  completed.  At  the  completion  of  the  connection  of 
the  two  models,  a  model  which  can  transform  finite-amplitude  waves  from 
deep  water  into  the  breaking  point  will  be  available. 

This  combined  cnoidal  II  -  Stokes  III  wave  model  will  be  best 
suited  for  open  coast  simulations,  modeling  monochromatic  waves  over 
smoothly  varying  bathymetries.  After  the  connection  research  is  com¬ 
pleted,  the  model  will  continue  to  be  improved.  The  comparison  of  model 
results  to  the  experimental  shoaling  data  indicate  that  a  third-order 
cnoidal  model  may  provide  an  improvement  over  second-order  cnoidal 
theory  near  the  breaking  point.  Development  of  the  necessary  third- 
order  cnoidal  wave  theory  expressions,  especially  for  the  energy  flux, 
would  be  an  arduous  task  very  well  suited  for  a  computer  program  using 
symbolic  manipulation. 

Additional  numer teal  techniques  and  wave  transformation  mechanisms 
will  be  added  to  expand  the  range  of  situations  for  which  the  model  is 
applicable.  The  incorporation  of  either  a  stretched  or  boundary  fitted 
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numerical  grid  would  enable  large  areas  to  be  more  economically  and 
accurately  modeled.  Perhaps  the  first  transformation  mechanism  which 
will  be  added  is  the  inclusion  of  energy  loss  due  to  bottom  friction. 

The  major  limitation  of  the  present  model  is  the  incapability  to  model 
highly  irregular  bathymetries  because  of  the  formation  of  caustics, 
regions  where  the  focusing  of  wave  energy  results  in  calculation  of 
unreliable  wave  properties.  Therefore,  the  addition  of  a  combined 
refraction-diffraction  calculation  scheme  would  be  a  significant 
improvement  that  would  greatly  add  to  the  versitility  of  the  model. 
Finally,  the  capability  to  model  a  spectrum  of  shallow  water  waves  must 
eventually  be  included.  However,  the  inclusion  of  a  spectral  capability 
will  require  extensive  research. 

Model  development  is  currently  hampered  by  the  lack  of  adequate 
laboratory  and  field  data  with  which  to  verify  numerical  wave  transfor¬ 
mation  models.  Excellent  shoaling  data  are  available  (and  used  in  this 
thesis),  but  adequate  data  for  refraction,  especially  data  for  wave 
angles,  are  not  available.  When  the  time  comes  to  verify  a  spectral 
shallow  water  finite-amplitude  wave  transformation  model  for  use  over  a 
complex  bathymetry,  if  adequate  data  are  not  available,  then  it  will  be 
difficult  to  design  the  optimum  model. 

The  ultimate  finite-amplitude  shallow  water  wave  transformation 
model  which  would  include  all  the  significant  wave  transformation 
mechanisms,  as  well  as  model  the  random  nature  of  the  sea  surface,  is 
not  clearly  in  sight.  It  is  hoped  that  this  report  may  be  a  step  in  the 
direction  of  that  final  goal. 
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APPENDIX  A 


EFFICIENT  CALCULATION  OF  ELLIPTIC  INTEGRALS 
AND  JACOBIAN  ELLIPTIC  FUNCTIONS 

This  appendix  contains  two  sections.  In  Section  A.1,  expressions 
for  the  efficient  calculation  of  elliptic  integrals  and  Jacobian  ellip¬ 
tic  functions  are  derived  based  on  the  work  of  Isobe  (1985).  Section 
A. 2  contains  a  short  table  of  properties  of  the  Jacobian  elliptic 
functions. 

A.1  Derivation  of  Expressions  for  the  Efficient  Calculation  of  Elliptic 
Integrals  and  the  Jacobian  Elliptic  Functions 

The  elliptic  modulus,  k  ,  is  a  fundamental  entity  which  enters  in 
the  definitions  of  the  elliptic  integrals  and  Jacobian  elliptic  func¬ 
tions.  As  is  discussed  in  Chapter  II,  <  is  used  as  the  auxiliary 
parameter  for  the  perturbation  expansion  leading  to  cnoidal  wave  theory. 
In  general,  ic  can  have  values  ranging  from  0  to  1  if  the  elliptic 
integrals  are  restricted  to  be  real-valued.  However,  in  cnoidal  water 
wave  theory,  ic  has  a  much  smaller  range.  At  U  =  25,  k  *  0.92  and,  as 
the  Ursell  number  increases,  k  approaches  1.0.  Computer  subroutines 
from  mathematics  libraries  and  formulae  found  in  references  such  as  Byrd 
and  Friedman  (1954)  and  Abramowitz  and  Stegun  (1972)  are  not  conven¬ 
iently  used  for  this  narrow  range  of  ic  values,  because  computer 
algorithms  based  upon  these  formulae  can  result  in  expensive  or  even 
inaccurate  calculations.  On  the  other  hand,  Isobe  (1985)  presented 
expressions  for  elliptic  integrals  and  Jacobian  elliptic  functions 
(Equations  3.6  to  3.12)  which  converge  quickly  and  accurately  for  the 
calculation  of  cnoidal  wave  theory.  Derivations  of  these  expressions 
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are  not  known  to  have  not  been  published  in  English.  Therefore,  they 


will  be  presented  in  this  appendix.  The  following  derivations  developed 


in  the  course  of  this  research  are  based  on  concepts  kindly  supplied  by 


Dr.  Isobe  (personal  communication,  1985). 


Several  definitions  are  necessary  for  the  derivations  that  fol¬ 


low.  The  complementary  modulus,  k'  ,  is  defined  from  the  following 


relationship: 


2  .2  i 

tC  +  tC  =1 


(  A.  1 ) 


The  complete  elliptic  integral  of  the  first  kind,  K  ,  and  its  comple¬ 


ment,  K'  ,  are  defined  as 


KU)  = 


_ _ 

2  2 

0  (1  -  k  sin2  *) 


(A. 2) 


K*  (k)  = 


(1  -  <'2  sin2*) 


(A. 3) 


The  complete  elliptic  integral  of  the  second  kind,  E  ,  and  its  comple¬ 


ment,  Ee  are  defined  as 


f  2  2  1/2 

E  (tc)  =  J  ( 1  -  k  sin  *)  d* 

e  0 


(A. 4) 


E’(k) 
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■/  < 


2  1/2 

k'  sin  *)  d* 


(A. 5) 
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The  nome,  q  ,  and  the  complementary  nome,  q'  ,  are  defined  as 


q 


-it(K'/K) 

e 


(A. 6) 


q' 


g-n ( K/K * ) 


(A. 7) 


Since  extensive  use  will  be  made  of  equations  from  Abramowitz  and 
Stegun,  (1972,  Chapters  16  and  17),  the  equation  numbers  from  this 
source  will  be  indicated  in  brackets  (eg.  [16.38.1])  for  easy 
reference.  Equations  necessary  for  the  calculation  of  elliptic 
integrals  and  the  Jacobian  elliptic  functions  are  readily  available. 
However,  these  expressions  are  often  given  in  terms  of  q  as  in 

(yj  '  =  1  ♦  2q  ♦  2q4  ♦  2q9  ♦  ...  [16.38.5]  (A. 8) 

=  2q1/4  (1  +  q2  +  q6  +  q12  +  ...)  [16.38.7]  (A. 9) 

(^T^)  =  (1  -  2q  +  2q4  -  2q9  ♦  ...)  [16.38.8]  (A.  10) 

Cnoidal  wave  theory  is  valid  for  U  5  10  (U  =  HL^/D^  ,  Equation 
2.23).  At  11  s  10  ,  q  and  q'  are  almost  equal  (q  »  q'  »  0.043). 

As  a  wave  moves  into  shallower  water,  the  Ursell  number  increases  and 
the  values  of  q'  and  q  become  smaller  and  larger,  respectively.  At 
U  =  100  ,  q  a  0.10  and  q'  a  0.00017  ;  at  U  =  500  ,  q  *  0.60  and 
q'  *  0.0000000039  •  Consequently,  for  the  range  of  Ursell  numbers 
applicable  to  cnoidal  theory,  expressions  using  power  series  based  on 
q'  will  converge  much  more  rapidly  than  expressions  based  on  q  .  The 

A3 


purpose  of  this  appendix  is  to  derive  expressions,  as  functions  of  q'  , 
which  are  needed  for  the  practical  and  efficient  calculation  of  cnoidal 
wave  theory. 

From  Equations  A. 2  and  A. 3,  it  can  be  seen  that 

K(k)  =  K'(k')  [16.1.2]  (A. 11) 

With  this  result  and  using  Equations  A. 6  and  A. 7,  Equations  A. 8  to  A. 10 
can  be  written  in  terms  of  q'  .  After  minor  algebraic  manipulations, 
the  results  are 


K’  =  |  (1  +  2q'  +  2q’4 


(A. 12) 


•c’  =  4q'1/2 
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(1  +  2q’  ♦  2q'4  ♦  ...) 


(A. 13) 


(A. 14) 


An  expression  relating  E  and  K  is  given  as 


-J  *  5  “  *  ',2)  *  (i)‘ 
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[17.3.23]  (A. 15) 


The  last  term  of  this  equation  can  be  simplified  using  a  binomial 


series.  The  result  is 
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(A. 16) 


As  was  the  case  for  K  and  K'  (Equation  A. 11),  E  and  E'  are 

C  v 

related  by 


E  (tc)  =  E*  (ic* ) 
e  e 


(A. 17) 


Substituting  Equation  A. 17  into  Equation  A. 16  and  using  Equation  A. 11 


yields 


e  1/i  2x 

jrr  =  ^  ( 1  +  <  )  + 


(k')  [l2 


2  (q'2  ♦  3q ' 


(A. 18) 


Equations  A. 12  and  A. 14  are  substituted  for  K'  and  k  on  the  right 
hand  side  of  Equation  A. 18.  After  carrying  out  the  multiplications, 
combining  terms,  and  neglecting  terms  of  q'^  and  higher  (in 
application  to  cnoidal  wave  theory,  the  largest  value  of  q'^  is 
approximately  0.0000005),  Equation  A. 18  becomes 


F '  2  4 

e  ( 1  +■  8q '  -  8q 1  ■*•...) 

K*  =  4  4 

(1  +  2q '  +  2q '  +  ...) 


(A. 19) 


Substituting  Equation  A. 12  into  the  left-hand  side  of  Equation  A. 19, 


E^  is  given  as 


*1(1  ♦  8q'2  -  8q 1 4) 


e  '  2 


(1  +  2q '  +  2q ’  ) 


(A. 20) 


umm 


Using  Equations  A. 7  and  A. 12,  an  expression  for  K  is 


v 1 

K  =  -  —  In  q' 

TT 

2 

=  -  (1  +  2*'  *  +  •••>  (A. 21) 


An  expression  for  Eg  in  terms  of  q  is 


E^K  +  KK'j 


[17.3.13]  (A. 22) 


1 

Substituting  expressions  for  Eg  ,  K'  ,  and  K  (Equations  A. 20,  A. 12, 
and  A. 21)  into  Equation  A. 22  results  in 


E 

e 


1  -  (-  j^-j(1  +  8q'2-  8q'4)  4-  ( 1  ♦  2q'  +  2q'4)4 


(1  +  2q'  +  2q’4) 


(A. 23) 


The  above  derivation  provided  expressions  in  terms  of  q'  for  the 
elliptic  integrals  (K,  K'  ,  Eg  ,  Eg),  the  elliptic  modulus  (k),  and 
the  complementary  modulus  (<')•  Next,  expressions  in  terms  of  q'  will 
be  derived  for  the  Jacobian  elliptic  functions;  cne  ,  sne  ,  dne  . 

Similar  to  the  relationship  between  the  cr.  and  the  cosine 
functions  (Chapter  III),  the  Jacobian  elliptic  functions  sn  and  dn  are 
analogous  to  the  sine  and  tangent  functions.  In  addition  to  the  three 
Jacobian  elliptic  functions  used  directly  to  specify  cnoidal  wave  theory 
( i . e . ,  cn,  sn,  and  dn),  three  other  Jacobian  elliptic  functions  are 
needed  in  the  derivations  which  follow.  These  quantities  can  be 
expressed  as  see  =  sne/cne,  nc9  =  1/cne  and  dc9  =  dn9/cn9  .  For  more 


information  the  reader  is  referred  to  Abramowitz  and  Stegun  (1972, 
Chapter  16). 

The  Jacobian  elliptic  functions  cn0,  sn0,  and  dn0  can  be  expressed 
in  terms  of  k'  ,  by  using  Jacobi's  imaginary  transform: 


sn(i0;tc)  =  iscCojK* ) 

[ 16.20. 1 ] 

(A. 24) 

cn(i0,<)  =  nc(Q;<' ) 

[16.20.2] 

(A. 25) 

dn(i0,<)  =  dc(0;ic' ) 

[  16.20.3] 

(A. 26) 

These  three  equations 

can  be  rewritten  using  the  following  relationship: 

0  =  i(-i0) 

(A. 27) 

Substituting  Equation 

A. 27  into  A. 24  to  A. 26  yields 

sn(0;ic)  =  isc(-i0,lc,) 

(A. 28) 

cn(0;«)  =  nc(-ie,x' ) 

(A. 29) 

dn(0;O  =  dc(-i0,<' ) 

(A. 30) 

A  general  expression  for  Jacobian  elliptic  functions 

is 

TO) 

pm.  = 

[16.36.3] 

(A. 31) 

In  which,  p  and  m 

can  be  any  of  the  symbols,  s  , 

c  ,  d  ,  or 

n  . 

The  three  elliptic  functions  of  interest,  cn0  ,  sn0  ,  and  dn  0  can  be 


fl 


expressed  in  terms  of  <'  using  Equations  A. 28  to  A. 31-  This  results 


sn(6;ic)  =  isc(-i9;«:' 


Hi 
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(A. 32) 


cn( 0 ;k  )  =  nc( -i0 ;tc ' ) 


Tn(-ie ' ) 

T  ( -  i  0 ;  < '  T 

c 


(A. 33) 


, 


dn( 0  ;< )  =  dc(-i0 ;< ' ) 


Td ( - i 9  ;  < ' ) 
TQ( -i0 ;k ' ) 


(A. 34) 


To  use  these  three  equations  to  evaluate  sn(0,q)  ,  cn(0,q)  and 

dn(0,q)  requires  expressions  in  terms  of  i0  and  q'  for  the  theta 
functions,  Tg  ,  TQ  ,  Td  ,  and  Tn  .  However,  these  expressions  are 
not  available  and  must  be  derived  from  the  following  equations  for  the 
theta  functions  given  in  terms  of  0  and  q  . 


Tg(0,q)  = 


iSr)  (si 


sin  S 


2  6  t 

-  q  sin  38  +  q  sin  58  + 


[16.38.1]  (A. 35) 


Tc(0;q)  = 
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cos  S 


2  6 

+  q  cos  36  +  q  cos  58  + 


'«>  =(fic)  ( 


1  +  2q  cos  28 


+  2q  cos  4e  + 


...) 


[16.38.2]  (A. 36) 


[16.38.3]  (A. 37) 


Tn(0;q)  = 


[2«'k) 


where 


1  -  2q  cos  28 
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[16. 38.4 ]  (A. 38) 


(A. 39) 


Using  the  following  relationships  between  trigonometric  functions  and 
hyberbolic  functions, 


cos  ( -  i 0 )  =  cosh  (-0)  =  cosh  8 


(A. 40) 


sin  ( - i 0 )  =  -  sin  i e  =  -  i  sinh  8 


(A. 41) 


Equations  A. 35  to  A. 38  can  be  expressed  in  terms  of  q'  .  This  results 
in  the  following  four  equations: 


Tg ( - i 9  ,* q ' )  =  -  if — 


^k^kK'  ^  S'  -  q'2  sinh  38'  +  q'^  sin  50'  +  . . 


(A. 1(2) 


/  1  /2  \  ^ 

r  (-i0;q')  =  f  2t'%, —  )  ( cosh  8'  +  q'2 


cosh  38'  +  q'  cosh  58' 


1  /  2 

Td( -ie ;q ’ )  =  +  2q'  cosh  2s'  +  2q cosh  116 '  +  "') 


J  \',d- 

(2K7)  T3^6' ^ 


(A. 44) 


/  \1/2  /  2i 

Tn(-ie;q')  =  (  ^K7/  v  ”  2q'  cosh  26'  +  2c*  cosh  ^6'  +  •  •  •) 


=  (2"K') 


T2(8') 


(A. 45) 


where 


,<  -  IL 


(A. 46) 


Tj(b')  =  sinh  0'  -  q'2  sinh  38'  +  q'^  sinh  58'  + 


(A. 47) 


r2(B')  =  1  -  2q'  cosh  28'  +  2q'  cosh  4e '  +  .. 


(A. 48) 


T^(8’)  =  1  +  2q'  cosh  28'  +  2q'  cosh  48'  +  .. 


(A. 49) 


2  6 

(s')  =  cosh  8'  +  q'  cosh  3s'  +  q'  cosh  58'  + 


(A. 50) 


The  next  three  expressions  will  be  used  to  simplify  the  derivation  which 


follows. 


Tq2  =  T2(0)  =  1  -  2q’  +  2q ' 


(A. 51) 


l  -  fcvtf 


T03  =  T  (0)  =  1  ♦  2q’  ♦  2q 


(A. 52) 


T4(0)  =  1  +  q'2  +  q'6  +  ... 


(A. 53) 


From  Equations  A. 32  to  A. 34,  A. 42  to  A. 45,  A. 13,  A. 14,  and  A. 51  to  A. 53, 
compact  notation  is  developed  for  cn(e,q)  ,  sn(8,q)  ,  and  dn(0,q). 

The  result  is 

sn(8,q)  =  isc  (-i0,q' ) 


=  i 


Tg  (-19, q’) 
Tc  (-19, q') 


(*) 


1/2  T1 (0* ) 


\  i2q'  +  2q'4\  T1(B’  } 
y1  -  2q '  +  2q'4/  T4(B' } 

T03  T1(fl,) 

T02  T4<6'> 


;  |9|  <  k 


cn(9,q)  =  nc(-i9,q’ ) 


Tn(-ie,q' ) 

Tc(-i9,q') 


*1/2  T2(B') 

v^t 


1  +q'2+q'6\T2(B,) 
J  -  2q'  +  2q'  y  T4(B' ) 

T04  T2(B,) 

T02  T4(0,) 


;  1 9 1  <  K 


A1 1 


(A. 54) 


(A. 55) 


dn(0,q)  =  dc(-ie,q' ) 


Td(-i9,q’) 

Tc(-ie,q’) 


/  «'  V/2  T3(S,) 

‘\4q'1/2y/  T4(S,) 

-7im|2m,6Vb(6,) 

\1  +  2q'  +  2q'/T4 


t3(b') 

VFT 


;  |6|  <  K 


(A. 56) 


With  these  expressions,  the  elliptic  functions  can  be  calculated  for 
| 0 |  >  K  using  the  following  three  expressions: 


cn(e  +  2nK;x)  =  (-1)ncn(9;ic)  (A. 57) 

sn(9  +  2nK;<)  =  (-1  )nsn(0,ic)  (A. 58) 

dn(0  +  2nK;ic)  =  dn(9;»c)  (A. 59) 

where 

n  =  0,1,2  — 


The  elliptic  integrals  can  be  written  using  this  compact  notation 
defined  in  Equations  A. 47  to  A. 53.  Thus,  Equations  A. 14,  A. 13,  A. 21, 
A. 12,  A. 23,  and  A. 20  become 


(A. 61) 


(A. 62) 


(A. 63) 

(A.6H) 

(A. 65) 

I 

(A. 66) 


These  expressions  (Equations  A. 54  to  A. 65)  are  the  equations  for  the 
efficient  calculation  of  cnoidal  wave  theory  first  presented  by  Isobe 
(1985).  These  equations  are  used  in  the  numerical  model  described  in 
this  thesis. 

A. 2  A  Table  of  Properties  of  Jacobian  Elliptic  Functions 

The  following  table  is  based  on  a  table  from  Isobe  and  Kraus  (1983b, 
Appendix  8).  It  is  included  here  for  easy  reference. 


Relations  between  squares  of  the  elliptic  functions: 


(A. 67) 


2  2 
sn  0  +  cn  e=1 


2  2  2 
k:  sn  0  +  dn  0  =  1 


(A. 68) 


Derivatives  of  the  elliptic  functions: 


-rr  (sn  0)  =  cn  0  dn  0 
d0 


—  (cn  6)  =  -  sn  0  dn  0 


d  .  2 

—  (dn  0)  =  -k  sn  0  cn  0 


Derivatives  of  even  powers  of  the  cn  function: 

_.2 


(A. 69) 


(A. 70) 


(A. 71) 


2n 


d0 


- 

l 

C 


(cn^“  0)  =  |2n(2n  -  1)xcn2n“2  0  +  (2n)c(1  -  xjcn^'1  0 


2n 


-2n(2n  +  1)cn2n+20 


) 


!^2X  + 


(cn2  0)  =  K2 | 2X  +  4(1  -  x)cn20  -  6cn^0 
d0 


] 


(cn^0)  =  k2  12xcn20  +  16(1  -  X)cn^0  -  2Ocn^0 

d0^  L 


(A.72) 


(A. 73) 


(A. 74) 


~  j7  (cn20)  =  ic4 1 8x ( 1  -  X)  +  8(2  -  1 3x  +  2X2)cn20 
d0 


-120(1  -  x)cn^0  +  12Ocn^0 J 


where 


X  = 


(1-0 


(A. 75) 


(A. 76) 
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The  derivations  of  the  equations  for  the  surface  profile,  n  , 
(Chapter  2);  the  energy  flux,  F  ,  and  energy,  E  ,  (Appendix  B)  require 
the  time-averaged  values  of  even  powers  of  the  cn  function.  Since  a 
cnoidal  wave  is  periodic  in  both  time  and  space,  then 


m. 
cn  9 


s  t/  °"m9  d9  ■  it 


m„  ,  „ 

cn  9  de 


(A. 77) 


From  Figure  3.2,  it  can  be  seen,  for  m  =  2,  4,  6...,  that 


it on"6  de  =  *t 


m .  . . 
cn  0  d0 


(A. 78) 


For  second-order  cnoidal  wave  theory,  the  calculation  of  three  of  these 
quantities  (for  m  =  2,  4,  and  6)  is  needed. 


Integrals  of  even  powers  of  the  cn  function: 


-K  _  rv/2  2  , 

cn2ede  =  f  - S23 _L 


f  cn2ed8  =  f 

Jc\ 


0  L  2  ,  2  A 

yl  -  k  sin  ^ 


d$ 


-££ 


-tf/2  .  2  .2 

1  -  <  sin  <fr 


d<j> 


2  .  2  , 
ic  sin  $ 


-l 


■tt/2 


1  -  K 


4 


d4> 


2  .  2  A 
1  -  <  sin  4> 


=  1-  E  -  1  K 


2  e 


(A. 79) 


From  Equation  A. 72, 


A 1 5 


f  cn2n+2  0  d0  =  |n-~  1  x  f  cn2n'2  0  d0 

I 0  2n  +  1  J0 


2n  f  1 


(1  -  X)  /  cn2n  0  d0 


(A. 80) 


Finally, 


—  2  1/2 

c,  =  cn  0  =  -  I  cn  0  d0  =  -X  +  y 

1  KJ0 


(A. 81) 


—  ii  i  u  i  ? 

=  cn  0  =  —  J  cn  0  d0  =  i  (-X  +  2y  +  2X  -  2Xy) 


(A. 82) 


—  6 

c,,  =  cn 


0  =  cn6  0  d0  =  -^  (-4x  +  8y  -  7Xu  +  3X2  +  8x2y  -  8x3) 


where 


(A. 83) 


(A. 84) 
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APPENDIX  B 


DERIVATION  OF  EXPRESSIONS  FOR  ENERGY  FLUX,  ENERGY,  AND 
GROUP  VELOCITY  FOR  SECOND-ORDER  CNOIDAL  WAVES 

This  appendix  contains  the  derivation  of  energy-related  expres¬ 
sions  needed  for  the  calculation  of  shoaling  and  refraction  of  second- 
order  cnoidal  waves.  Corresponding  results  for  first-order  cnoidal 
waves  are  derived  in  the  process.  Section  B.1  contains  a  derivation  of 
an  expression  for  energy  flux  of  second-order  cnoidal  waves.  Section 
B.2  contains  a  derivation  of  an  expression  for  the  average  energy  per 
unit  surface  area  of  second-order  cnoidal  waves.  Finally,  in  Section 
B.3  the  derivations  of  the  first  two  sections  are  used  to  derive  an 
expression  for  the  group  velocity  of  second-order  cnoidal  waves. 

B.1  Energy  Flux  of  Second-Order  Cnoidal  Waves 

From  Phillips  (1977,  p.  63),  the  average  energy  flux,  F,  in  the 
direction  of  wave  propagation,  per  unit  surface  area,  is 

F  =  p  u  ^  (u2  +  w2)  +  gz  +  ^  dz  (B.1) 

The  overbar  denotes  averaging  over  wavelength  as  in  Appendix  A  (Equation 
A. 77).  Note  that  energy  flux  is  a  vector  quantity,  but  since  F  is 
defined  in  the  direction  of  wave  propagation,  the  vector  notation  has 
been  omitted  for  convenience.  Direct  evaluation  of  this  equation  using 
second-order  cnoidal  wave  expressions  for  the  variables  would  be  tedi¬ 
ous.  However,  a  less  complicated  expression  for  F  can  be  derived  (Isobe 
and  Horikawa,  1982).  Using  the  Bernoulli  equation  for  unsteady  motion, 


(B.2) 


pi  ,2  2,  x  PB 

“  +  5  (u  +  w  )  +  gz  -  4>.  =  — 

P  d  t  p 


Equation  B.1  becomes 


dz 


(B.3) 


For  waves  of  permanent  form,  the  horizontal  wave  particle  orbital 
speed,  u,  can  be  expressed  as 


u  =  u(x  -  Ct,z) 


(B.i*) 


Since  u  =  and  w  =  ,  the  most  general  expression  for  the 


velocity  potential,  <j>  ,  is 


4>  =  4»(x  -  Ct,z)  +  f(t) 
=  $(X,z)  +  f(t) 

where  X  =  x  -  Ct  . 

Evaluating  the  time  derivitive  of  $  ,  using 


^t  =  *XXt  +  ft 


-c*x  +  ft 


and 


^  -  $ „X 

x  TX  x 


it  can  be  seen  that 


>t  =  Cu  +  ft 


Sustituting  Equation  B.8  into  Equation  B.3  yields 


F  =  p 


/  u(Cu*r*r t) 


dz 


B2 


(B.5) 


(B.6) 


(B.7) 


(B.8) 


(B.9) 


It  can  be  shown  that  f(t)  can  be  chosen  such  that  the  Bernoulli 
constant,  p^/p  i  vanishes  and  such  that  f^.  is  not  a  function  of  t 
(see,  eg.,  Isobe  and  Kraus,  1983b,  Appendix  C).  Therefore,  since 
Stokes'  second  definition  of  wave  celerity  (Equation  2.21)  was  used  to 
derive  the  second-order  cnoidal  wave  expressions  (Chapter  II),  Equa¬ 
tion  B.9  simplifies  to 


F  = 


( B. 10) 


For  use  in  the  numerical  model,  this  equation  must  be  expressed  in 

terms  of  the  second-order  cnoidal  wave  theory  presented  in  Chapter  II. 

The  power  series  expansion  of  F,  based  on  the  perturbation  parameter, 

e  =  H/D  (Equation  2.44),  will  be  developed.  The  expression  for  F  will 

include  the  first  two  terms  of  the  perturbation  power  series.  In  the 

course  of  the  derivation,  these  two  approximations  will,  at  first, 

2  3 

appear  to  be  0(e  )  and  0(e).  However,  the  final  form  of  the  equation 

for  F  will  show  these  first  two  approximations  to  be  0(1)  and  0(e).  The 

2 

third  contribution  to  the  power  series  of  F  (coefficients  of  e  )  is  not 
included  because  the  third  approximation  is  not  complete  for  a  second- 
order  theory.  This  point  will  be  illustrated  during  the  derivation. 
Repeating  Equation  2.67,  u  can  be  expressed  as 

u  =  [B00  *  B10°t  *  B20°2 

p  4 

-sH"2)  <B01  *  B1t°1  *  B2,C2>J  (B',1) 


The  terms  in  this  equation  must  be  expanded  and  grouped  into  factors  of  a 
common  power  of  e  .  From  Equations  2.75  to  2.80,  Equation  B.11  becomes 


u  =v4b,  *  6^b2  "  ('  *  ‘S)C1  *  'Vs 

-iM2( 

=  |^e ( b1  +  c1 )  +  e2(fc>2 

- !  M 


2  2  2 
e  b5  +  e  b6C1  +  e  b 


7°s)] 


*  °i>  *  '  (b2  *  b3=,  *  b4c2> 


e  (b5  +  b6c1  +  b 


7C2)  ] 


where 


b1  =  X  -  u 


b2  =  ^  (X  -  y  -  2X2  +  2y2) 


b3  =  jj  ( 1  -  6X  +  2y ) 


b4  =  -1 


h  -  3x 
b5  =  "2 


b6  =  3  -  3X 


7  ‘  2 


(B. 12) 


(B. 13) 


(B. 14) 


(B. 15) 


(B. 16) 


(B. 17) 


(B. 18) 


(B. 19) 


Before  making  substitutions  to  write  Equation  B.12  in  a  more  com¬ 


pact  form,  it  will  improve  the  clarity  of  the  following  manipulations  to 


examine  the  variables  in  Equation  B.12.  The  only  z  dependency  is  con¬ 


tained  in  the  (z+D)/D  term.  The  only  terms  which  vary  with  wavelength 


are  c1  and  c2-  The  variables  H  and  D  are  assumed  not  to  vary  over 


wavelength  as  a  basic  assumption  in  the  derivation  of  waves  of  permanent 
form  in  Chapter  II.  The  b's,  functions  of  elliptic  functions  (and, 
therefore,  the  Ursell  parameter),  are  constant  over  both  depth  and 
wavelength. 

Rewriting  Equation  B.12  in  a  more  compact  form  results  in 

u  =  yfi D  +  £2  (u2  -  \  Y2  U3)]  (B.20) 

where 

U,  =  b1  +  c1  (B.21 ) 

U2  —  b2  +  b^c  ^  ^  ^||^2  ( B .  22 ) 

U3  =  b^  +  b^c^  +  b?c2  (B.23) 

Y  =  2-±-2  (B.24) 

Squaring  u  and  gathering  coefficients  of  like  powers  of  e  yields 

u2  =  gD  £e2U2  +  e3(2U1U2  -  Y2) 

♦  e4  (u2  -  U2U3  Y2  +  £  U2  Y4)]  (B.25) 


Referring  back  to  Equation  B.10,  uc  must  now  be  integrated  over  depth  and 
then  averaged  over  wavelength.  The  resulting  equation 


B5 


can  be  simplified  to 


2  .  n  /  2  ,2.  3  ,  2S  4  .  2X 

u  dz  =  gD  l  e  <u^>  +  e  <u^>  +  e  <u^> 

The  three  integrals  on  the  right-hand  side  of  Equation  B.26  are 
~2  ~2  ~2 

represented  by  <u^>  ,  <u^>  ,  and  <u^>  ,  respectively.  The  symbols,  <  > 
are  used  to  denote  integration  over  depth,  from  -D  to  n. 

From  Equation  B.10,  it  is  seen  that  Equation  B.27  must  be 
multiplied  by  pC  .  The  expression  for  wave  celerity,  C,  repeated  from 
Equation  2.65,  is 


(B.27) 


c  =  Vir(Co  +  eC1  ♦  e2C2)  (B.28) 

Now,  using  Equations  B.27  and  B.28,  B.10  can  be  written 

nFy»  S  e2C  <u2>  +  e3(c  <U2>  +  C  <u2>\ 

PgD  -ygD  o  1  \  O  2  1  1  / 

♦  e4  ^Cq  <u2>  C1  <u2>  +  C2  <u2>^  +  HOT  (B.29) 

As  was  previously  mentioned,  only  the  first  two  terms  of  Equation  B.29 

are  complete  for  a  second-order  theory.  This  can  be  seen  by  studying 

the  equation  for  u  (Equation  B.12).  A  third-order  theory  would  con- 

3  ? 

tribute  terms  to  this  equation  which  are  coefficients  of  e  .  For  u  , 

H 

the  lowest  terms  unique  to  the  third-order  theory  would  be  at  the  e 

•3  i  h 

level  (e  x  e  )  .  Therefore,  the  coefficient  of  e  in  Equation  B.29 
is  not  complete  for  a  second-order  theory.  Keeping  only  the  complete 
orders,  Equation  B.29  is  rewritten  as 


B6 


lW 


(B.30) 
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l 

s 
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i 

l 

■fc 


pgDT/gTT =  £"  co  <ui>  +  e3(co  <U2>  +  c i  <ui>)+  hot 


Although  it  appears  that  only  two  orders  remain  in  Equation  B.30 

2  3  — 2 

[coefficients  of  0(e  )  and  0(e  )],  this  is  not  the  case.  Both  <u^> 

— p 

and  <u^>  are  expressed  as  power  series  of  e  .  Therefore,  to  maintain 

only  two  levels  of  approximation  in  the  evaluation  of  Equation  B.30,  it 

~2  ~2 

is  necessary  to  truncate  the  contributions  of  <u^>  and  <u2>  .  For  the 
2 

coefficient  of  e  in  Equation  B.30,  only  the  lowest  two  orders,  [0(1) 

~ 2  3 

and  0(e)],  of  <u^>  will  be  allowed.  For  the  coefficient  of  e  only  the 

— p  ~2 

lowest  order,  [0(1)],  of  both  <u^>  and  <u2>  will  be  allowed. 

~5 

From  Equations  B.27,  B.26,  and  B.21  <u.>  can  be  written  as 


-r 


(b1  +  c ^ )  dz 


2  2 
^  +  2b^c^  +  c^  dz 


(B.31) 


Since  there  is  no  z  dependency  in  this  equation,  the  depth  integration 

~ 2 

is  easily  accomplished.  The  result  of  the  depth  integration  of  <u^>  is 


2S  ( |2 

nbi 


+  2b1c1  +  c^  j  (n  +  D) 


(B.32) 


An  expression  for  n  is  needed  to  evaluate  this  equation.  From  Equation 


2.66,  n  can  be  written  as 


n  =  D  ^e(a1 


+  c1)  *  e  (a2  +  a3C1  +  ¥2) 


(B. 33) 


/s-.‘ss‘s- 


where 


(B.34) 


a1  =  X  -  y 


2  '  4 


.  +  y  -  2X^  +  2Xy) 

(B. 35) 

t  .  .  3 

3  "  4 

(B.36) 

a  -  1 
a4  '  4 

(B.37) 

The  only  terms  which  vary  over  wavelength  are  and  c2;  noting  from 

Equation  A. 81  to  A. 83  that  c2  =  and  c^  =  =  c^,  then 

Equation  B.32  becomes  1 


'1 


=  °[' 


<u?>  =  D  I  b?  +  2b1c1  +  c^  +  e  (a^  +  2a^b^c 


1 


ai°t  *  biet  *  2bi°2 


?  *  4 


+  HOT  ( B . 38 ) 


2  2 
Evaluating  <u2>  is  more  complicated  than  the  evaluation  of  <up  because 

~5 

<u2>  has  Y-terms,  which  are  functions  of  depth.  From  Equations  B.27, 
B.26,  and  B.21  to  B.23,  <u2>  can  be  written  as 


<u2>  = 


■  f  h 


+  c1}  (b2  +  b3°1  +  W 


-  (b,  +  c,)  (bg  +  bgC,  +  b?c2)  _dz 


( B . 39 ) 


Carrying  out  the  depth  integration  yields 


<u2>  =  2( b 1  +  c 1 )  (b2  +  b^c1  +  b4c2^  +  D> 


-  (b1  ♦  c  1 )  (b5  *  bgC^  +  b?c2)  (n  +  ~)+  H0T 


(B.40) 


B8 


This  equation  can  be  simplified  by  neglecting  higher-order  contributions. 
~2 

Since  <U2>  is  included  only  in  the  term  representing  the  second 

~2 

approximation  of  F  (Equation  B.30),  only  terms  of  0(1)  in  <U2>  will 
survive  in  the  final  expression  for  F.  Inspecting  the  equation  for 
n  (Equation  B.33),  it  can  be  seen  that  multiplication  by  n  in  Equation 
B.40  will  increase  the  order  of  the  result.  The  terms  resulting  from 
the  multiplication  of  n  are  of  higher  order  and  are  dropped,  and 
Equation  B.40  becomes 

<u^>  =  2  (b1  +  c^ )  (b2  +  b3c1  +  b^)  D 

-  (b1  +  c 1 )  (b  +  bfic^  +  b?c2)  |  (B.41 ' 

carrying  out  the  multiplication  this  equation  becomes 

<u2>  ■  2D(b1b2  *  b,b3V  *  blVl  *  b2V  *  b3°1  *  b4c1  ) 

-  I(b1b5  *  blVT  *  b1b7°?  *  Vi  *  b6c1  *  b7c1  )  (B-,2> 

Finally,  all  of  the  terms  in  the  expression  for  F  (Equation  B.30)  are 
ready  to  be  assembled.  Substituting  Equations  B.38  and  B.42  into  B.30 
gives 


pgD  ^gD* 


coDV 

f  L  2 

Vb1  + 

*  «3[c( 

»D  (a1 

.  2DC0 

(  b1b2 

-  5  c 

3  0 

(b1b5 

>i°i  *  A ) 


2  —  2  2  2  3  \ 

>1  +  2a^b^c^  ♦  a^c^  +  b^c^  +  2b^  +  ci  ' 

.  2DC„  (  b.b,  .  b,b  ^  .  b,b,c2  .  b^  .  b  e2  .  b(C3  ) 


2  .  — 


1  6  1  17  1  5  1  6  1 

♦  C. 


>7°?) 


:,D(b2  ♦  2b, C,  .  O2  ) 


(B.43) 
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For  use  in  the  numerical  model,  this  expression  for  F  is  put  in  terms 
of  X  and  y  .  Substituting  Equations  B.34,  B . 1 3  to  B . 1 9  ,  A. 81  to  A. 83, 
and  2.65  into  Equation  B.43  yields 

- pF  —  =  e2  (-X  ♦  2y  +  4Xy  -  X2  -  3w2)  4 

PgD^VgD  J 

+  e3  ( — 4X  +  8y  +  53Xu  -  12X2  -  60y2 

+  53X2y  -  120Xy2  -  8x3  +  75u3)  ^  (B.44) 

The  apparent  order  of  this  equation  can  be  reduced  by  using  the  following 
relationship  based  on  the  definition  of  e  (Equation  2.44) 

H2  =  D2e2  (B.45) 

Substituting  this  equation  into  Equation  B.44  results  in  the  final 
expression  for  F.  The  energy  flux  for  first-order  cnoidal  wave  theory 
is  proportional  to  FQ;  whereas,  the  eF  term  is  the  contribution  for 
second-order  cnoidal  wave  theory. 

F  =  H2  ogyfgD  (FQ  ♦  zF}) 

=  H2  og^gD"  ( -X  -*■  2y  +  4Xy  -  X2  -  3y2)  ^ 

♦  e  (-  4x  ♦  8y  ♦  53Xu  -  12x2  -  60y2 

53x2y  -  120Xu2  -  8X3  ♦  75u3)  ^ 


♦ 


(B.46) 


B.2  Energy  of  Second-Order  Cnoidal  Waves 


The  average  wave  energy  per  unit  surface  area,  E,  can  be  expressed  as 
the  sum  of  the  average  potential  energy,  Ep  ,  and  the  average  kinetic 
energy,  Ek  . 


E  =  Ek  +  Ep  (B.47) 

The  two  contributions  in  this  equation  will  be  evaluated  in  the  next  two 
sections. 


B.2.1  Kinetic  Energy 


The  equation  for  average  kinetic  energy  per  unit  surface  area  is 


u 


dz  + 


( B . 48 ) 


The  quantities  necessary  for  expressing  the  first  integral  in  this 
equation  were  developed  in  the  previous  section.  From  Equation  B.27, 
the  depth  integration  of  u2  averaged  over  wavelength  is 


(B.49) 


As  is  the  case  for  the  derivation  of  F,  E  will  be  expressed  up  to  and 

~2 

including  the  second  approximation.  Equation  B.30  shows  that  <u^>  will 

~2 

contain  terms  of  0(1)  and  0(e)  and  <u^>  will  contain  terms  of  0(1). 


Assembling  these  contributions  from  Equation  B.38  and  B.42  and  inserting 
them  into  Equation  B.49  results  in 

/*n  o  il  t  i  _  o  v 


£  u2  dz  =  gD2L2  (b2  ♦  2b, ^  *  c2  ) 

*  '3  [a1b?  *  2a1b1=i  *  a1e?  *  1 


c,  *  a,c2  *  b2c,  *  2b, c2  *  0 ] 


+  2 

(b1b2  * 

b1b3°1 

+  b\bnc*  + 

b2C1 

+  b^c^  +  b^c3 

1 

'  3 

.  .  ~2 

,  ~2  ,  ~3 

(b1b5  * 

blVT 

+  b1b?c1  + 

Vi 

+  b6c1  +  Vi 

-  3  (b1b5  *  blVl  *  blb7°1  *  Vi  *  Vi  *  VWJj  (B'50) 

Expressing  this  equation  in  terms  of  X  and  ti  and  using  Equation  B.45 
yields 


P  u2  dz  =  gH2  ^  (-X  +  2u  *  4Xy  -  X2  -  3n2) 

+  e  ( X  -  2y  -  2Xy  +  3X2 
-  15p2  -  2X2y  -  30Xy2  +  2X3  +  30y3 )  ~ 


(B.5D 


The  second  integral  of  Equation  B.48  will  require  an  evaluation  similar 
to  the  evaluation  of  u2  given  in  the  derivation  of  F.  Repeating 
Equation  2.68,  the  vertical  water  particle  orbital  speed,  w,  is  given  as 


"  ■  ViD«coad  (2_S_S)  ( B10  *  2B20°l) 

-^(H-^Vll  *2B21bl)  < 

From  Equations  2.76  to  2.80,  B.13  to  B . 1 9 ,  and  B.24,  Equation  B.52 
becomes 

w  =  y[gD  4K  £  csd  j  eY  +  e2  Jy^  +  2bljc1) 


(B.52) 


5  y3  (b6  +  2b 


7C1 '] 


(B.53) 


B12 


This  equation  written  in  more  compact  form  is 


w  =  yfgD  4K  £  csd  ^eY  +  e2  (  YV^  -  g-  Y3W2  )J 


where 


W1  =  b3  +  2Vl 


W2  +  b6  +  2b7C1 


(B.54) 


(B.55) 


(B.56) 


Squaring  w  and  gathering  coefficients  of  like  powers  of  e  results  in 

w2  =  gD  16K2  (csd)2  [\2Y2  +  e32y(yW1  -  £  Y3W?  ) 

♦  *"  (  ™,  -  5  2 )  ] 


(B.57) 


Equation  B.57  must  be  integrated  over  depth  and  then  averaged  over 
wavelength.  The  resulting  equation 


/■ 


2  dz  =  gD  16K2  £  (csd)2 


e2  /  Y2  dz 


/ 


/n 

21  (™,  -  5“2)dz  *  '4  /  (™i 

-D 


can  be  simplified  to 


J r  w2  dz  =  gD  16K2  (csd)2  ^  e2  <w‘ 


I  -2)  * 

(B.58) 


3  .  2.  4.2. 

+  e  <w2>  +  e  <w^> 


(B.59) 


where  the  three  integrals  on  the  right-hand  side  of  Equation  B.58  are 

~2  ~2  ~2 
represented  by  <w1>,  <w2>,  and  .w^>. 

~5  ~2 

Just  as  <u^>  contributed  only  higher-order  terms  to  F,  <w^> 
contributes  only  higher-order  terms  to  E^.  Futhermore,  because  of  the 

p  ~2 

term  (D/L)  in  Equation  B.59,  it  will  be  shown  that  only  <w^>  will 
contain  terms  of  sufficiently  low  order  to  contribute  to  E^.  From 
Equations  2.23  and  2.44  (U  4  e),  the  following  equation  relates  the 
Ursell  parameter  and  the  perturbation  parameter. 


(c/  ■  5 


(B . 60) 


Since  the  smallest  value  of  U  for  which  cnoidal  wave  theory  is  valid  is 
approximately  10,  (D/L)2  <  e.  This  implies  that  (D/L)2  contains  a 
"hidden"  e,  which  will  increase  the  order  of  Equation  B.59.  Substituting 
Equation  2.85  (dispersion  relation)  into  Equation  B.60  results  in 


.  _3 
'  16 


2„2 

K  K 


e  +  0(  e  ) 


(B.61) 


This  expression  is  substituted  into  Equation  B.59  resulting  in  Equation 
B.62. 


{ 


w2  dz  =  gD  l6K2y|  ~2~2  (csd)2  <w^> 
K  K 


4  2  5  2 

♦  e  <wp  +  e"*  <w*> 


(B.62) 


Note  that  Equation  B.62  is  one  order  higher  than  Equation  B.27,  the 


I 

) 

corresponding  equation  for  the  u2  component  of  Ek.  Maintaining  a 
consistent  contribution  from  both  components  requires  that  only  the 
lowest  order  of  Equation  B.62  remain.  Thus  Equation  B.62  becomes 


{ 


w2  dz  =  (csd)2  e3  <w2> 

K 


(B.63) 


Therefore,  only  the  lowest  order  term  of  <w*>  will  be  represented. 

2 

Referring  to  Equations  B.24,  B.59  and  B.58,  <w^>  can  be  expressed  as 


■r 


<w2>  =  /  Y2  dz 


=  ^  +  HOT 


(B.64) 


Substituting  this  equation  into  Equation  B.63  yields 


w2  dz  =  (csd)2  e3 

tc 

Using  Equation  B.45,  Equation  B.65  can  be  written  as 


n  _ 

w2  dz  =  H2  (csd)2  e 
< 

Recall  (from  Equation  2.84)  that  the  definition  for  csd  is 


(B.65) 


(B.66) 


csd  =  cn  0  sn  0  dn  0  (B.67) 


where,  cn,  sn,  and  dn  are  Jacobian  elliptic  functions  and  0  is  defined 


B15 


in  Equation  4.2.  These  functions  do  not  depend  on  z,  but  they  are  a 
function  of  x.  Therefore,  they  must  be  averaged  over  wavelength.  After 
squaring  Equation  B.67,  the  result  can  be  expressed  in  quantities  for 
which  the  average  over  wavelength  is  known.  Using  the  following  two 
identities  repeated  here  from  Equations  A. 67  and  A. 68. 

sn2e  +  cn2e  =  1  (B.68) 

K2sn20  +  dn2e  =1  (B.69) 


the  square  of  Equation  B.67  can  be  expressed  as 


2  2  22  P  U  ii  P  f\ 

(csd)  =  cn  0  -  <  cn  0  +  2ic  cnH0  -  cn  o  -  k  cna0 


=  C1  -  +  2k%  *  c2  ’  KS 


(B.70) 


Substituting  this  equation  into  Equation  B.66  yields 


/' 


2  2  /  °1  —  —  C2  — 

w  dz  =  gH  |  -  c1  +  2c„ - p  -  c.  |  e 

K  K  * 


( B .  7 1 ) 


The  definition  of  X  (rewritten  from  Equation  2.70) 


x  = 


1  -  K 


(B.72) 


can  be  used  to  simplify  the  terms  of  the  above  equation  which  have  a  <c 


in  the  denominator 


(B.73) 


11 

2 


-  Xc^  +  C,| 


2  =  Xc2  +  c2 


(B.74) 


resulting  in 


w2  dz  =  gH2  (  xc 


y  +  c2  -  Xc2  -  c^  ^  e 


(B.75) 


This  equation  is  put  in  terms  of  X  and  y  using  Equations  A. 81  to  A. 83. 


/n 

w2  dz  =  gH2  (-X  +  2y  +  2Xy  -  3X2  +  2X2y  -  2X3j  e  (B.76) 

Assembling  the  two  components  of  Ek  (Equations  B . 5 1  and  B.76)  results  in 
Equation  B.77,  the  average  kinetic  energy  per  unit  surface  area  of  a 
second-order  cnoidal  wave. 


dz 


\-2  /  “2  d*  *  / 

=  pgH2  h  (-X  +2 p  +  4Xu  -  X2  -  3p2) 

+  e  ^-X  +  2Xy  +  2y  -  3X2  -  15y2 


+  2X2y  -  30  Xu2  -  2X3  +  30y3 


)ss] 


B.2.2  Potential  Energy  of  Second-Order  Cnoidal  Waves 


The  average  potential  energy  per  unit  surface  area  is 


(B.77) 


#  r ’i 

Ep  =  pg  z  dz 


(B.78) 
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The  depth  integration  of  this  equation  results  in  two  terms. 

E*  =  \  pgn2  +  \  PgD2  (B.79) 

The  second  term  in  this  equation  represents  the  potential  energy  of  the 
still  water  column.  Since  we  are  interested  in  the  portion  of  the 
potential  energy  due  to  the  waves,  only  the  first  term  of  the  above 
equation  will  be  evaluated.  Therefore,  the  potential  energy  per  unit 
surface  area  due  to  waves  is 


E 

P 


(B.80) 


Squaring  n  (Equation  B.33)  and  keeping  only  the  first  two  approximations 
results  in 


n2  =  D2  [(a2  ♦  2alCl  ♦  c2  )  e2  ♦  2  (a^  ♦  aia3Cl 

+  a^c2  +  a2c1  +  a^c2  +  a^c3  je3J+  HOT  (B.81) 

Making  use  of  Equations  B.45,  B.34  to  B.37,  and  A. 81  to  A. 83,  Equation 
B.80  becomes 


Ep  =  PgH 


[i( 


•X  +  2p  ♦  4Xy  -  X2  -  3v»2 


)*  e(x  • 


2p 


1 2 X y  +  3x2  +  5u2  -  12X2u  +  10Xy2  +  2X3 


(B.82) 


B.2.3  Expression  for  Energy  of  Second-Order  Cnoidal  Waves 


Combining  Ek  (Equation  B.77)  and  Ep  (Equation  B.82)  results  in  Equation 
B.83,  the  equation  for  the  average  wave  energy  per  unit  surface  area. 


As  for  F,  the  contribution  for  first-order  wave  cnoidal  theory  is 


proportional  to  Eq,  and  the  contribution  for  second-order  cnoidal  wave 
theory  is  proportional  to  eE^. 


E  =  pgn  (E0  +  eE1 ) 


x  +  2u  +  4xy  -  x2  -  3v2) 


=  pgH2  [l(-x  +  2y  +  4xy  -  x2  -  3y2) 

♦  e  (x  -  2y  -  17Xp  +  3X2  -  17X2y  +  2X3  +  15y3) 


(B.83) 


B.3  Group  Velocity  of  Second-Order  Cnoidal  Waves 


The  velocity  of  wave  energy  propagation;  the  group  velocity,  C  ,  is 

© 

defined  (Phillips,  1977,  p.  25)  as 


C  =  ^ 
g  '  E 


(B.84) 


Using  the  results  for  F  and  E  developed  above  (Equations  B.46  and  B.83), 


Equation  B.84  becomes 


pgH2ygD*  (Fq  +  eF^ ) 

g  pgH2  (Eq  +  eE1 ) 


(B.85) 


The  following  algebraic  manipulations  on  Equation  B.85  will  lead  to  a 
simplified  form  for  Cg. 


/Fn  +  eFi 

Ce 


.  WFo  ^  eFi 


E°  \  1  El 

v  +  e^ 


WFo ;  £F1 


Eo  1  Ei 

v  +  e  *S 


/  E1  \  2  E1 

E°  ,  /E'  M  2/EA2 

L  ,**U**s  r‘U)  . 


.  &  *  o(t2> 

Eo  1 


(B.86) 


From  Equations  B.H6  and  B.83, 


E0  =  F0 


(B.87) 


The  final  form  for  Cg  is  obtained  by  substituting  Equation  B.87  into 
B.86  yielding 


c*  =  VIM1 


Fi  -  Ei 


(B.88) 


Note  from  the  expression  for  C,  Equation  B.28,  that  C  and  Cg,  at  the 
first  approximation,  are  both  equal  to  the  shallow  water  wave  speed  of 
small -amplitude  wave  theory,  Vi^*- 


vyi'/r, 


I  V  *  '  » 
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COMPUTER  PROGRAM 


PROGRAM  CNOIDAL 


PARAMETER  <NX=110 ,NY  =  1 10 ,NZ=NX*NY ) 

REAL  LAM,MU,L,M 

COMMON  /Al/  D(NY,NX>,  H(NY,NX>,  LAM( NY ) ,  MU (NY) 

COMMON  / A2 /  L ( NY , NX ) ,  C'NY.NX),  C0(NY>,  Cl ( NY ) ,  C2(NY), 

5  ALFA( NY ,NX) 

COMMON  /A 3/  Q ( NY , NX  > ,  U(NY,NX>,  DELQMAXf  NX ) 

COMMON  / A4 /  F  (  NY ,  NX ! ,  FO(NY),  FI i NY ) ,  Eone ( NY ) ,  Cg(NY.NX) 
COMMON  /A5/  II, NI, JJ.NJ.IT, III 
COMMON  / A6 /  P I , DELHMAX , GAM , G , T , DXDY , A 


CHARACTER  TITLE*20,  ORDER* II,  READDPT*3,  READBC*3 ,  BEACH*5 , 
$  RED0F*3 

DIMENSION  ITER(NX) ,LBP(NY> 


C  INNTIALIZE  VARIABLES 

DATA  Q/NZ*0. 000000001- 
DATA  PI  /  3. 141592654/ 


C  READ  INPUT  FILE 

OPEN  (UNIT=1 ,NAME=  INPUT.DAT  , STATUS =' OLD ' ) 

READt 1,801)  A, ALPHA, HO .T.G.RHO 

READ( 1,801)  DX.DY.SU.TOL 

READ( 1,803'  II ,NI , JJ ,NJ . I OUT, III 

READ* 1,004)  READDPT. BEACH, D1 ,B,M 

READ! 1,805)  READBC , NT 


DXDY  *  DXDY 

GAM  =  RHO+G 

ALPHA  *  ALPHA*PI ' 180. 


READ  DEPTHS  OR  CALCULATE  DEPTHS 
ALT  DATUM  CORRECTION  'OR  SURGE) 

IF- READDPT< 1:3) .EQ.  YES  )  THEN 


DO  5  1=1. NI 

READ. 1,806)  ■ D( J.I > ,J=1,NJ' 

DO  5  J=1 ,NJ 

D-  J . I  )  =  D- J  ,  I  ♦  SU 

CONTINUE 


PNY  *  FLOAT-  NJ  -  JJ 
PNX  =  FLOAT  NI 

c i :  =  pi  *  : .  - 


DO  -  I  =  1  .  N I 
PI  ■  FLOAT  I 
1ELX  =  PI  *  DX 
XX  -  PI  RNX 


do  6  t=j:.nj 

t  *  =  FL " AT  ’  TT 
YY  =  P ’  PNY 

C"  •  P  »  M*DELX 

"NTINUE 


Dl*XX*COS '  PI2*YY 


D1*XX 


ENI  IF 


C  SET  BOUNDARY  CONDITIONS  ON  SEAWARD  BOUNDARY. 

IF(READBC< 1:3) .EQ.  YES' )  THEN 
READ( 1,955)  <H( J,NI) ,J*JJ,NJ) 

READ! 1 , 960 )  (ALFA! J,NI) ,J=JJ,NJ> 

DO  20  J=JJ,NJ 

ALFA! J,NI)  *  ALFA! J,NI )  *  PI/180.  +  PI 
20  CONTINUE 

ELSE 

DO  30  J*1 ,NJ 
ALFA!J,NI>  *  ALPHA 
H! J,NI )  =  HO 
30  CONTINUE 

END  IF 

DO  12  J=JJ ,NJ 
H!J,NI-1)  *  H( J ,NI ) 

12  CONTINUE 

C  CALCULATE  ELLIPTIC  QUANTITIES  ON  SEAWARD  BOUNDARY. 

MM  =  5 

CALL  ELLIP(NI.MM) 

C  CALCULATE  CELERITY  AND  WAVE  LENGTH  ON  SEAWARD  BOUNDARY. 
CALL  LENGTH (NI) 

C  CALCULATE  ENERGY  FLUX  ON  SEAWARD  BOUNDARY. 

NF  =  0 

CALL  EFLUX(NI,NF) 

NF  =  1 

C  CALCULATE  TRUE  URSELL  NUMBER  ON  BOUNDARY 
DO  50  J=JJ ,NJ 

U(J,NI)  =  H( J,NI )  *  L( J ,NI ) **2  /  D( J ,NI ) **3 
50  CONTINUE 

C  START  OF  LOOP  TOWARD  SHORE  FROM  BOUNDARY  OR  CONNECT  POINT. 
C  LOOP  ON  ROWS  OF  CONSTANT  I. 

DO  300  I=NI , II+l , -1 
IM  =  1-1 

C  ITERATION 

DO  200  IT=1 , 20 
ITER(IM)  =  IT 
DELHMAX  =0.0 
DELQMAX(IM)  =0.0 

C  CALCULATE  ELLIPTIC  QUANTITIES. 

MM  =  2 

CALL  ELLIP(IM,MM) 

C  CALCULATE  CELERITY  AND  WAVE  LENGTH. 

CALL  LENGTH (IM) 

C  SOLVE  FOR  WAVE  ANGLE 
CALL  ANGLE! IM) 

C  CALCULATE  WAVE  HEIGHT 

CALL  EFLUXt IM,NF) 

C  DO  LATERAL  BOUNDARIES 

ALFA( JJ , IM)  =  ALFA(JJ+1,IM) 

ALFA(NJ,IM)  =  ALFA(NJ-l.IM) 

F( JJ.IM)  =  F<  JJ+1 , IM) 

F(NJ, IM)  =  F(NJ-1,IM) 

HfJJ.IM)  =  H( JJ+1 , IM) 

H(NJ.IM)  =  H(NJ-1 , IM) 

IF( DELHMAX.LT. TOD  GOTO  201 

200  CONTINUE 

201  CONTINUE 


C  2 


C  “BREAKING"  CRITERION 
C  PLANE  BEACH  BREAKING  CRITERION 

IF< BEACH' 1:5>.EQ. 'PLANE' >  THEN 

HD  *  H< ♦ , IM)  /  D( 4 , IM) 

IF< HD. GE. BREAK)  THEN 
III  »  IM 

HDLAST  •  H( 4 , I )  /  D'4,I> 

HDDIFF  ■  HD  -  HDLAST 

HD1  «  (HD  -  BREAK)  /  HDDIFF 

HD2  «  'BREAK  -  HDLAST)  /  HDDIFF 

H(4, IM)  *  HD1*H< 4.1)  ♦  HD2*H'4.IM> 

D( 4 , IM)  *  HD1*D<  4,1)  ♦  HD2*D'4,IH> 

ALFA( 4 , IM)  «  HD1*ALFA< 4,1)  ♦  HD2*ALFA< 4 , IM) 

XB  «  FLOAT ( IM )  *  OX  ♦  HD1  *  DX 

GOTO  301 

ENDIF 

ELSE 

C  IRREGULAR  BOTTOM 

HSTOP  ■  10 

DO  250  J ■JJ  ,NJ 

HD  -  H(J.IM)  /  D<  J , IM) 

IF( HD. GT. BREAK)  THEN 

IF( IIII. EQ. 10000)  IIII  -  IM 

H(J.IM)  •  D( J , IM)  *  BREAK 

IF(LBP< J) .EQ.O)  LBP(J)  «  IM 

REDOF  -  YES 

ENDIF 

HD  *  H'J.IM)  /  D'J.IM) 

HSTOP  -  AMINK HSTOP. HD) 

250  CONTINUE 

IF ( HSTOP. GE. BREAK)  THEN 
III  -  IM 
GO  TO  301 
ENDIF 

C  IF  HEIGHTS  ARE  TRUNCATED  BY  BREAKING  CRITERION  THEN  RECALCULATE  F. 
IF( REDOF' 1:3) .EQ.  YES  >  THEN 
CALL  ELLIP< IM.MM) 

CALL  LENGTH* IM) 

NF-0 

CALL  EFLUXI IM.NF) 

NF-1 

ENDIF 

ENDIF 

C  SET  HAVE  HEIGHT  ON  NEXT  I  ROW  TO  HEIGHT  AT  CURRENT  ROW  IN 
C  PREPARATION  FOR  NEXT  ITERATION 
C  ALSO  CALCULATE  TRUE  'JR SELL  NUMBER  FOR  OUTPUT. 

DM  •  IM- 1 

DO  275  J-JJ.NJ 
HIJ.IW)  -  H'J.IM) 

V'J.IM)  »  H(J.IM)  *  LiJ.IM>**2  i  D<  J , IM) **3 
275  CONTINUE 

300  CONTINUE 

301  CONTINUE 


c 


OUTPUT 


DO  325  I-II.NI 
DO  325  J-JJ.NJ 

ALFA(J.I)  *  ( ALFA < J , I )  -  PI)*180./PI 
325  CONTINUE 

ALPHA  *  (ALPHA  -  PI)*180./PI 

C  SET  HAVE  HEIGHTS  INSIDE  BREAKING  POINT  TO  ZERO  FOR  OUTPUT 
DO  400  I  *  l.NI 
DO  400  J  *  JJ  ,NJ 

IF(LBP(J) .GT. I)  THEN 
H( J , I )  «  0. 

ALFA ( J ,  I )  *  0. 

L( J.I)  *  0. 

Cg  <  J , I )  -  0. 

U(J.I)  -  0 . 

END  IF 

400  CONTINUE 

IF(A.EQ.O.O)  THEN 
ORDER  *  ' CNOIDAL  I' 

ELSE  IF(A.EQ.l.O)  THEN 
ORDER  -  CNOIDAL  II' 

ELSE 

ORDER  -  BONKERS 
ENDIF 

WRITE* 11,810)  ORDER 
HRITE( 11,815) 

HRITE( 11.820)  HO. T. ALPHA 

WRITE* 11.825)  OX.DY.TOL 

HRITE( 11.830)  SU.G.READBC 

WRITE' 11.832)  READDPT.B.M.Dl , III ,NI 

WRITE* 11.815) 

WRITE* 11.835)  *  ITER* I > , I - 1 1 ,NI ) 

I®  ITE ( 11*815) 

WRITE* 1 1  *  840 )  *  DELQMAX < I > , I ■ I I , NI ) 

WITE*  11.815) 

TITLE-  WAVE  HEIGHTS 
CALL  POUT* 100.0. H, TITLE) 

TITLE-  HAVE  ANGLES 

CALL  POUT* 10.0, ALFA. TITLE) 

IF*  IOUT.GE.  1)  THSf 

TITLE-  DEPTHS 

CALL  POUT* 100. ,D. TITLE) 

ENDIF 

IF* IOUT.GE. 2)  THEN 
TITLE-  URSELL  NUMBER 
CALL  POUT* 1.0. U, TITLE) 

ENDIF 

IF* IOUT.GE. 3)  THEN 
TITLE-  HAVE  LENGTH 
CALL  POUT* 10. ,L. TITLE) 

ENDIF 


IF( IOUT.GE. 4 )  THEN 
TITLE* ' HAVE  CELERITY' 
CALL  POUTUO.  ,C, TITLE) 
ENDIF 


IF( IOUT.GT. 5)  THEN 
TITLE* ' ENERGY  FLUX' 

CALL  POUT(0. 01, F, TITLE) 
ENDIF 

IF( IOUT.GE. 6)  THEN 
TITLE* ' GROUP  VELOCITY' 
CALL  POUT< 100. ,Cg, TITLE) 
ENDIF 


801  FORMAT (6F1 0.1) 

803  FORMAT(6I10) 

804  FORMAT (A6 ,T11 ,A5,T21 , 3F10 . 5) 

805  FORMAT( A6 , Til , 110 ) 

806  FORMAT (8F1 0.4) 

810  FORMAT( / ,T20, All, 'REFRACTION  AND  SHOALING') 

815  FORMAT ( / ) 

820  FORMAT ( '  INPUT  WAVE  *  ' ,F6. 4, T30, 'WAVE  PERIOD  =  ',F6.3,T55, 
$  '  INPUT  ANGLE  =  ' .F6.1) 


825  FORMAT ( '  DELTA  X  * 

$  '  TOLERANCE  = 

830  FORMAT ( '  SURGE 
$  '  READ  BC 

832  FORMAT ( '  READ  HEIGHTS 
$  '  SLOPE  = 

$  '  AMPLITUDE  * 

$  '  PRINT  STOP  * 


'  ,F6.2,T30, 'DELTA  Y  = 

' ,F6.5> 

'  ,F6.2,T30,'G 
'  ,A3 ) 

A3, T30,’ INTERCEPT  = 

' , F6 .4 ,  / , 

'  ,F6. 3, T30,  PRINT  START  * 
'  ,16,/) 


'  ,F6. 2,T55, 


'  ,F6.3,T55, 

' ,F6.3,T55, 


'  ,I6,T55, 


835  FORMAT( 4013 ) 

840  FORMAT (10G1 0.3) 

905  FORMAT( 15X,8F12.8) 

910  FORMAT ( 315 ,8F12 .8 ,F12 . 5 ) 
950  F0RMAT(A8) 

955  FORMATt  10F8. 3 ) 

960  F0RMATU0F8.1) 

END 


) 

l 


C  5 


c  ********************************************************************* 

SUBROUTINE  POUT ( FACT , DUM , TITLE ,  IOUT  > 

PARAMETER  (NX-110, NY-110, NZ-NY*NX> 

COfttoN  /A5/  II,NI,JJ,NJ,IT,III 
DIMENSION  KX(NY) ,DUM(NY,NX) 

CHARACTER  TITLE* 16 ,LINE*5 

LINE-' . ’ 

WRITE( 11,900)  TITLE 
HRITE< 11,905)  FACT 

NC  -  20 

JMAX  -  NJ  -  JJ  +  1 
JT  »  JMAX/NC 
JTT  -  MOD ( JMAX, NC) 

IF(JTT.NE.O)  JT-JT+1 
J1  =  JJ 
J2  -  Jl+NC-1 
IF( J2.GT.NJ)  J2-NJ 

DO  100  KK-l.JT 

HRITE< 11,910)  (J,J-J1,J2) 

HRITE( 11,912)  (LINE, J-Jl ,J2) 

DO  50  I-III ,NI 
DO  40  J-J1.J2 
RND  »  0.5 

RND  -  SIGN< RND,DUM< J ,1) ) 

KX(J)  -  INT ( FACT*DUM( J , I )  +  RND) 

40  CONTINUE 

HRITE( 11,915)  I , (KX( J) , J-Jl , J2) 

50  CONTINUE 
C  WRITE( 11 ,920) 

J1  »  J1  +  NC 
J2  -  J2  +  NC 
IF< J2.GT.NJ)  J2-NJ 

100  CONTINUE 

900  FORMATS, T2.A16) 

905  FORMAT ( T2, ' (MULTIPLIED  BY  ',F10.2,')') 

910  FORMAT( / , '  I/J  s ' ,20( 2X,I3) ) 

912  FORMAT ( 6X , 2  0A5  > 

915  FORMAT( IX, 13, IX, ' i ' ,2015) 

920  FORMAT < ' 1 ' ) 

RETURN 

END 


P  *AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
SUBROUTINE  ELLIP(I.MM) 

PARAMETER  <NX=110 ,NY=110 ) 

REAL  LAM, MU 

COMMON  /Al/  D(NY,NX) ,  H(NY,NX) ,  LAM (NY) ,  MU(NY) 

COMMON  /A3/  Q(NY,NX) ,  U(NY,NX) ,  DELQMAX(NX) 

COMMON  /A5/  II,NI,JJ,NJ,IT,III 
COMMON  /A6/  PI,DELHMAX,GAM,G.T,DXDY,A 

C  STATEMENT  FUNCTIONS  FOR  THE  CALCULATION  OF  ELLIPTIC  QUANTIES 

S(J,I)  =  1.0  +  8.0*Q< J,I)**2  -  8 . 0*Q( J, I ) **4 

T02( J, I)  =  1.0  -  2 . 0*Q( J , I )  +  2.0*Q( J,I)**4 

T03( J,I)  =  1.0  +  2.0*Q< J,I)  +  2 . 0*Q( J , I ) **4 

T04( J,I)  *  1.0  +  Q(J,I)**2  +  Q( J,I)**6 

DO  101  J*JJ,NJ 

el  *  H(J,I)/D( J.I) 

U(J.I)  =  G*H( J,I)*T*T  /  <D( J,I)*D( J.I) ) 

C  USE  SUCCESSIVE  UNDERRELAXATION  PARAMETER  (EXPONENTIAL  FIT). 

R  *  EXP  <  -0.61452  -  0.05249*U(J,I)  ) 

DO  30  M-l.MM 

LAM(J)  *  16 . 0*Q  ( J ,  I )  *  ( T04 (J.I)  /  T02(J,I>)**4 
Q1  ■  A*  <  -1.0  -  2.0  *  LAM(J)  )  /  4.0 

QQ  .EXP  (  -SQRT  (  0. 75*U< J, I ) *( 1 .  ♦  el*Ql>  '  >  T02(J,I>**2  ) 
DEL Q  *  00  -  Q(J.I) 

OOLD  -  Q(J.I) 

0(J,I)  -  QQ  -  (DELQ*R) 

30  CONTINUE 

DELQPC  -  al>s<  (  QOLD  -  Q(J,I)  )  )  /  QOLD 
DELQMAX( I >  »  AMAX1 ( DELQMAX ( I ) , DELQPC ) 

C  CALCULATE  NECESSARY  ELLIPTIC  QUANTITIES. 

LAM(J)  -  16. 0*Q( J, I )*(T04<  J, I )/T02(J,I) )*M 

MU(  J)  «  (  2.0/  (-ALOG  (  Q(  J,  I )  )  )  -  S(J,I)  ♦  T03(J,I)*M)  / 

$  T02( J,I)**4 

101  CONTINUE 

RETURN 

END 


■  'j' 


Avr  i"™  V.Y -n«  V, *»**•* ;»y»Vir 

;  .'■  *  1  /  y'  .  •  ■  ~ 1 


0 

SUBROUTINE  LENGTH  (I) 

PARAMETER  ( NX-110 ,NY*110 ) 

REAL  LAM, MU, L 

COMMON  /Al/  D(NY,NX> ,  H(NY,NX),  LAM(NY),  MU(NY) 

COMMON  /A2/  L(NY,NX),  C*NY,NX),  CO(NY),  Cl ( NY) ,  C2(NY), 

$  ALFA (NY, NX) 

COMMON  /A5/  II ,NI , JJ,NJ , IT, III 
COMMON  / A6/  PI , DELHMAX . GAM , G , T , DXDY , A 

DO  102  J *JJ,NJ 

el  *  H(J,I)/D<J,I) 
e2  *  el*el 

CO ( J )  *  SQRT*G*D< J.I) ) 

Cl ( J '  *  (  1.0  ♦  2 . 0*LAM( J )  -  3 . 0*MU < J )  )*0.5 
C2<J>  -  A  *  (  -  6.0  -  16 . Q*LAM( J )  +  5.0*MU(J> 

$  -  16 . 0ALAM<  J )  **2  10 . 0*LAM(  J  )  *MU<  J  ) 

9  ♦  15 . 0*MU ( J  S  **2  )  t  40.0 

C*J,I)  *  C0*J)*<  1.0  ♦  el*Cl < J )  «■  e2*C2( J)  ) 

L ( J ,  I )  *  C< J.I)*T 

102  CONTINUE 

RETURN 

END 

C  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 
SUBROUTINE  ANGLE* I ) 

PARAMETER  ( NX-110 .NY* 110 ) 

REAL  L 

comou  IK2I  L*NY,NX>.  C(NY.NX),  CO(NY),  C1*NY>,  C2(NY), 

9  ALFA* NY, NX) 

COMMON  /A5/  II, NI.JJ.NJ.IT, III 
COMMON  /A6/  PI. DELHMAX, GAM. G.T. DXDY, A 

IP  -  I+l 

DO  45  J* JJ+-1 . NJ - 1 
JM  -  J-l 
JP  -  J*1 

ALFA *  J ,  I  >  -  ASIN(  <SIN<  ALFA(J.IP)  >  *  L(J.IP)  -  0.5*DXDY* 

9  '  COS *  ALFA*  JP . IP ) )  '  L(JP.IP)  -  COS ( ALFA* JM, IP > )  '  L(JM.IP)  )) 

(  *  L'J.I)  ) 

C  CORRECT  FOR  AS IN  FUNCTION 

ALFA* J, I)  •  PI  -  ALFA* J , I ) 

45  CONTINUE 

RETURN 

END 


lii. 


V? 


C8 


c  ********************************************************************* 

SUBROUTINE  EFLUX(I,NF) 

PARAMETER  <NX*110, NY-110) 

REAL  LAM, MU, L 

COMMON  /Al/  D(NY,NX) ,  H<  NY,NX) ,  LAM (NY) ,  MU(NY) 

COf«ON  /A2/  L(NY,NX),  C(NY,NX) ,  CO(NY),  Cl ( NY) ,  C2(NY), 

$  ALFA (NY, NX) 

COMMON  /A4/  F(NY,NX) ,  FO(NY) ,  Fl(NY),  Eone(NY),  Cg(NY,NX> 

COttlON  /AS/  II ,NI , JJ ,NJ, IT, III 
COMMON  /A6/  PI ,DELHMAX,GAM,G,T,DXDY,A 

DO  103  J*JJ,NJ 
RL  »  LAM(J) 

RM  *  MU( J ) 

RL2  *  RL*RL 

RL3  *  RL2*RL 

RM2  *  RM*RM 

RM3  *  RM2*RM 

el  »  H(J,I)  /  D( J , I ) 

F0( J)  *  (  -RL  +  2.0*RM  ♦  4.0*RL*RM  -  RL2  -  3.0*RM2  )  /  3.0 
F1(J)  *  A  *  (  -4 . 0*RL  +  8 . 0*RM  +  53.0*RL*RM  -  12.0*RL2 

9  -60 . 0*RM2  +  53 . 0*RL2*RM  -  120.0*RL*RM2 

9  -8. 0*RL3  *  75.0*RM3  >  /  30.0 

Eone(J)  *  A  *  (  RL  -  2.0*RM  -  17.0*RL*RM  +  3.0*RL2 
9  -17 . 0*RL2*RM  ♦  2.0*RL3  +  15.0*RM3  )  /  30.0 

Cg2  *  (  F1(J)  -  Eone(J)  )  /  F0(J) 

Cg(J.I)  *  C0(J)  *  (  1.0  ♦  el  *  A  *  Cg2  ) 

103  CONTINUE 

IF(NF.EQ.O)  THEN 
DO  104  J  *JJ ,NJ 
el  *  H(J,I)  /  D< J,I) 

F3  *  GAM  *  C0(J)  *  (  F0(J)  *  el*Fl(J)  ) 

F( J ,  I )  -  H( J , I ) **2  *  F3 

104  CONTINUE 
ELSE 

IP  «  I+l 

DO  105  J-JJ+l.NJ-1 
JP  -  J+l 
JM  »  J-l 

el  »  H< J,I)  /  D(J,I) 

F3  -  GAM  *  CO<J)  *  (  F0< J)  ♦  el*Fl(J>  ) 

C  SOLVE  FOR  ENGERY  FLUX  (MX  DOT  F) 

F( J , I )  •  (  F(J.IP)  *  COS ( ALFA ( J , IP) )  ♦  0.5  *  DXDY  * 

9  <  F ( JP , IP )  *  SIN(ALFA( JP,IP) )  -  F(JM,IP)  *  SIN( ALFA( JM, IP) )  )  ) 

9  /  COS(ALFA(J,I) ) 

HNEH  «  SQRT<  F(J.I)/F3  ) 

DELPC  •  ABS(  HNEH  -  H(J,I)  )  /  H(J,I) 

H( J ,  I )  «  HNEH 

OELHMAX  -  AMAXKDELHMAX, DELPC) 

105  CONTINUE 
END  IF 

RETURN 

END 


C  9 


APPENDIX  D 


REFRACTION  OF  SECOND-ORDER  CNOIDAL  WAVES  OVER  A  PLANE  BOTTOM  -  PLOTS 

This  appendix  contains  plots  of  nondimension  wave  height  vs.  depth 
for  the  Refraction  of  second-order  cnoidal  waves  over  a  plane  bottom 
slope  of  1:50.  The  plots  contained  in  this  appendix  are  companion  plots 
to  Figures  6.15  to  6.29,  which  are  plots  of  nondimension  wave  angle  vs. 
depth.  For  complete  information  on  these  simulations  the  reader  is 
referred  to  Chapter  VI,  Section  6.2.2. 
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Figure  D.1  Wave  Height  -  Refraction  over  a  plane 
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Figure  D.2  Wave  Height  -  Refraction  over  a  plane 
bottom  (1:50),  a  =  30°,  T  =  8  s 
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Figure  D.4  Wave  Height  -  Refraction  over  a  plane 
bottom  (1:50),  aQ  =  30°,  T  =  12  s 
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APPENDIX  E 


NOTATION 

The  following  symbols  are  used  in  this  thesis: 

a  Depth  at  shoreline  for  calculated  bathymetry 

an  Term  used  to  simplify  n  (n  =  0,1, 2, 3, 4) 

Ajj  Term  in  equation  for  n  (n  =  0,1,2) 

b  Bottom  slope  of  plane  beach  portion  of  calculated  bathymetry 

bn  Term  used  to  simplify  u  and  w  (n  =  1,2, 3... 7) 

b00  Coefficient  in  zeroth-order  of  the  stream  function 

Bmn  Term  in  equations  for  u  and  w  (m  =  0,1,2  ;  n  =  0,1) 

cn  cnne  (n  =  1,2,3) 

C  Wave  celerity 

Cg  Magnitude  of  group  velocity 

C  Group  velocity 

© 

Cn  Nth-order  component  of  wave  celerity  (n  =  0,1,2...) 

dn  Nth-order  component  in  expansion  of  (D/L)  (n  =  1,2,3...) 

D  Water  depth  measured  to  the  Stillwater  level 
d  Amplitude  of  cosine  portion  of  rhythmic  bathymetry 

E  Average  wave  energy  per  unit  surface  area 

En  Nth-order  component  to  E  (n  =  0,1) 

Ee  Complete  elliptic  integral  of  second  kind 
E;  E(k') 

E^  Average  kinetic  energy  per  unit  surface  area  due  to  waves 
Ep  Average  potential  energy  per  unit  surface  area  due  to  waves 

E*  Ep  plus  potential  energy  of  still  water 


f(t) 

Function  added  to  Bernoulli  equation  to  make  Bernoulli 
constant  vanish 

F 

Average  wave  energy  flux  in  direction  of  wave  propagation 

"F 

Average  wave  energy  flux 

'  f 

Fn 

Nth-order  contributions  to  F  (n  =  0,1) 

FS4I 

Wave  energy  flux  from  Sarpkaya  and  Isaacson  (1981) 

f3,i 

Component  of  energy  flux 

8 

Acceleration  of  gravity 

H 

Wave  height 

,v 

V 

Hb 

Wave  height  of  breaking 

Ho 

Deepwater  wave  height 

i 

=  \ pr 

*V 

a 

i 

Unit  vector  in  x-direction 

>  ■ » 

3 

Unit  vector  in  y-direction 

k 

Wave  number,  2n/L 

kSA 

Wave  number  from  small-amplitude  wave  theory 

'  • 

K 

Complete  elliptic  integral  of  the  first  kind 

K' 

K(k') 

'‘,<1 

Kr 

Refraction  coefficient  for  small -amplitude  waves 

r*.<v 

Ka 

Shoaling  coefficient  for  small -amplitude  waves 

L 

Wavelength 

'  , 

lsa 

Wavelength  from  small-amplitude  wave  theory 

1  *■ 

Lo 

Deep  water  wavelength 

— 

m 

Number  of  nodes  in  x-direction 

/-  * , 

V,’ 

..  «.  •< 

n 

Number  of  nodes  in  y-direction 

N 

Nondimens ional  water  surface  elevation 

Nn 

Nth-order  contribution  to  N  (n  =  1,2,3...) 

— 

E2 

_ 

p  Pressure 

pB  Bernoulli  constant 

P  Nondimens ional  pressure 

Pg  Nondimens ional  Bernoulli  constant 

Pn  Nth-order  component  to  expansion  of  P  (n  =  0,1,2...) 

q  Nome  of  the  theta  function  or 

Flow  rate  (used  only  in  Chapter  II) 

q'  Complementary  nome  of  the  theta  function 

* 

q'  Value  of  q  before  relaxation 

Q  Nondimens ional  flow  rate 

Qn  Nth-order  component  to  expansion  of  Q  (n  =  0,1,2...) 

R  Under-relaxation  parameter 

S  Quantity  used  in  calculation  of  E 

t  Time 

T  Have  period 

T^( )  Neville's  notation  for  theta  functions  (i  =  c,d,n,s) 

Tj^O  Theta  functions  (i  =  1,2, 3, 4) 

T0i  T(0)  (i  =  2,3,4) 

u  Horizontal  component  of  the  orbital  water  particle  velocity 
u  Oribital  water  particle  velocity 

p 

<un>  Term  used  to  simplify  calculation  of  E  and  F  (n  =  1,2,3) 

U  Ursell  parameter 

Un  Term  used  to  simplify  calculation  of  F  (n  =  1,2,3) 

w  Vertical  component  of  the  orbital  water  particle  velocity 
<wn>  Terra  used  to  simplify  calculation  of  E  (n  =  1,2,3) 

Wn  Term  used  to  simplify  calculation  of  E  (n  =  1,2,3) 

x  Horizontal  coordinate 
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xb 

Location  of  wave  breaking 

Ax 

Node  spacing  in  x-direction 

X 

Nondimens ional  horizontal  coordinate  or 
=  x-Ct  (in  Equations  B.5  to  B.1  only) 

Ay 

Node  spacing  in  y-direction 

1 

y 

horizontal  coordinate 

1 

Y 

=  (z  +  D)/D 

1 

z 

Vertical  coordinate 

Z 

Nondimensional  vertical  coordinate 

'v 

a 

Angle  of  wave  propagation 

.  f 

,  t' 

°b 

Wave  angle  at  breaking 

a 

o 

Deepwater  wave  angle 

•M* 

6 

=  ir8/2K 

■  *  t* 

y  %  \ 

•  ’  * 

8' 

=  ir0/2K' 

i 

6 

Auxiliary  parameter  for  the  perturbation  expansion 

»«;*, 

$ 

e 

=  H/D,  perturbation  parameter 

W. 

n 

Water  surface  elevation 

nSA 

n  as  determined  from  small  amplitude  wave  theory 

,  * 

nn 

Nth-order  components  of  n  (n  =  1,2,3-.-) 

e 

Argument  to  the  Jacobian  elliptic  functions 

0SA 

Small -amplitude  wave  phase  function 

l  » 

K 

Modulus  of  the  elliptic  integral 

•'/  i 

K' 

Complementary  modulus  of  the  elliptic  integral 

. 

X 

Term  used  to  simplify  elliptic  calculations 

u 

Term  used  to  simplify  elliptic  calculations 

t  t 

c 

The  elliptic  integral 

It 

The  constant  w 

— 

EH 

n 

p 

♦ 

« 

t 

T 

n 

Q 

V 

<> 


The  unified  nonlinearity  parameter  of  Goda 
Density  of  water 

Dummy  argument  of  elliptic  integral  or 
Velocity  potential 

Nondimens ional  velocity  potential 

Stream  function 

Nondimens ional  stream  function 

Nth-order  contribution  to  the  stream  function  (n  =  0,1,2...) 

The  portion  of  the  wave  phase  used  in  deriving  the  wave  angle 
equation 

Horizontal  gradient  operator 

Denotes  integration  over  depth  from  -D  to  n 

Denotes  integration  over  wavelength 
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